Francois Brun / gammalib

Commit 5fa82e7999a730f9ef0668ce7bf1fcee5dbdd9b0

Authored by knodlseder
1 parent 5f9063bf

update

Inline Side-by-side
Showing 16 changed files with 250 additions and 142 deletions
include/GEvent.hpp
  View file @5fa82e7
... ... @@ -20,6 +20,9 @@
20 20 #define GEVENT_HPP
21 21  
22 22 /* __ Includes ___________________________________________________________ */
  23 +#include <string>
  24 +#include <iostream>
  25 +#include "GLog.hpp"
23 26 #include "GInstDir.hpp"
24 27 #include "GEnergy.hpp"
25 28 #include "GTime.hpp"
... ... @@ -60,6 +63,10 @@
60 63 ***************************************************************************/
61 64 class GEvent {
62 65  
  66 + // I/O friends
  67 + friend std::ostream& operator<< (std::ostream& os, const GEvent& event);
  68 + friend GLog& operator<< (GLog& log, const GEvent& event);
  69 +
63 70 public:
64 71 // Constructors and destructors
65 72 GEvent(void);
... ... @@ -69,19 +76,18 @@ public:
69 76 // Operators
70 77 virtual GEvent& operator= (const GEvent& event);
71 78  
72   - // Event access methods
  79 + // Pure virtual methods
  80 + virtual void clear(void) = 0;
  81 + virtual GEvent* clone(void) const = 0;
  82 + virtual double size(void) const = 0;
73 83 virtual const GInstDir& dir(void) const = 0;
74 84 virtual const GEnergy& energy(void) const = 0;
75 85 virtual const GTime& time(void) const = 0;
76 86 virtual double counts(void) const = 0;
77 87 virtual double error(void) const = 0;
78   -
79   - // Other methods
80   - virtual void clear(void) = 0;
81   - virtual GEvent* clone(void) const = 0;
82   - virtual double size(void) const = 0;
83   - virtual bool isatom(void) const = 0;
84   - virtual bool isbin(void) const = 0;
  88 + virtual bool isatom(void) const = 0;
  89 + virtual bool isbin(void) const = 0;
  90 + virtual std::string print(void) const = 0;
85 91  
86 92 protected:
87 93 // Protected methods
... ...
include/GEventAtom.hpp
  View file @5fa82e7
... ... @@ -20,6 +20,7 @@
20 20 #define GEVENTATOM_HPP
21 21  
22 22 /* __ Includes ___________________________________________________________ */
  23 +#include <string>
23 24 #include "GEvent.hpp"
24 25 #include "GInstDir.hpp"
25 26 #include "GEnergy.hpp"
... ... @@ -64,19 +65,20 @@ public:
64 65 // Operators
65 66 virtual GEventAtom& operator= (const GEventAtom& atom);
66 67  
67   - // Event access methods
  68 + // Pure virtual methods
  69 + virtual void clear(void) = 0;
  70 + virtual GEvent* clone(void) const = 0;
  71 + virtual double size(void) const { return 1.0; }
68 72 virtual const GInstDir& dir(void) const = 0;
69 73 virtual const GEnergy& energy(void) const = 0;
70 74 virtual const GTime& time(void) const = 0;
71   - double counts(void) const { return 1.0; }
72   - double error(void) const { return 0.0; }
  75 + virtual double counts(void) const { return 1.0; }
  76 + virtual double error(void) const { return 0.0; }
  77 + virtual std::string print(void) const = 0;
73 78  
74 79 // Other methods
75   - virtual void clear(void) = 0;
76   - virtual double size(void) const { return 1.0; }
77   - virtual GEventAtom* clone(void) const = 0;
78   - virtual bool isatom(void) const { return true; }
79   - virtual bool isbin(void) const { return false; }
  80 + bool isatom(void) const { return true; }
  81 + bool isbin(void) const { return false; }
80 82  
81 83 protected:
82 84 // Protected methods
... ...
include/GEventBin.hpp
  View file @5fa82e7
... ... @@ -20,6 +20,7 @@
20 20 #define GEVENTBIN_HPP
21 21  
22 22 /* __ Includes ___________________________________________________________ */
  23 +#include <string>
23 24 #include "GEvent.hpp"
24 25 #include "GInstDir.hpp"
25 26 #include "GEnergy.hpp"
... ... @@ -66,19 +67,20 @@ public:
66 67 // Operators
67 68 virtual GEventBin& operator= (const GEventBin& bin);
68 69  
69   - // Event access methods
  70 + // Pure virtual methods
  71 + virtual void clear(void) = 0;
  72 + virtual GEvent* clone(void) const = 0;
  73 + virtual double size(void) const = 0;
70 74 virtual const GInstDir& dir(void) const = 0;
71 75 virtual const GEnergy& energy(void) const = 0;
72 76 virtual const GTime& time(void) const = 0;
73 77 virtual double counts(void) const = 0;
74   - virtual double error(void) const;
  78 + virtual double error(void) const = 0;
  79 + virtual std::string print(void) const = 0;
75 80  
76 81 // Other methods
77   - virtual void clear(void) = 0;
78   - virtual double size(void) const = 0;
79   - virtual GEventBin* clone(void) const = 0;
80   - bool isatom(void) const { return false; }
81   - bool isbin(void) const { return true; }
  82 + bool isatom(void) const { return false; }
  83 + bool isbin(void) const { return true; }
82 84  
83 85 protected:
84 86 // Protected methods
... ...
include/GammaLib.hpp
  View file @5fa82e7
... ... @@ -124,9 +124,14 @@
124 124 #include "GModelPar.hpp"
125 125 #include "GModelSpatial.hpp"
126 126 #include "GModelSpatialPtsrc.hpp"
  127 +#include "GModelSpatialConst.hpp"
  128 +#include "GModelSpatialCube.hpp"
127 129 #include "GModelSpectral.hpp"
128 130 #include "GModelSpectralPlaw.hpp"
  131 +#include "GModelSpectralPlaw2.hpp"
129 132 #include "GModelSpectralExpPlaw.hpp"
  133 +#include "GModelSpectralFunc.hpp"
  134 +#include "GModelSpectralConst.hpp"
130 135 #include "GModelTemporal.hpp"
131 136 #include "GModelTemporalConst.hpp"
132 137  
... ...
include/Makefile.am
  View file @5fa82e7
... ... @@ -29,9 +29,14 @@ pkginclude_HEADERS = GammaLib.hpp \
29 29 GModelPar.hpp \
30 30 GModelSpatial.hpp \
31 31 GModelSpatialPtsrc.hpp \
  32 + GModelSpatialConst.hpp \
  33 + GModelSpatialCube.hpp \
32 34 GModelSpectral.hpp \
33 35 GModelSpectralPlaw.hpp \
  36 + GModelSpectralPlaw2.hpp \
34 37 GModelSpectralExpPlaw.hpp \
  38 + GModelSpectralFunc.hpp \
  39 + GModelSpectralConst.hpp \
35 40 GModelTemporal.hpp \
36 41 GModelTemporalConst.hpp \
37 42 GApplication.hpp \
... ...
inst/cta/include/GCTAInstDir.hpp
  View file @5fa82e7
... ... @@ -20,7 +20,7 @@
20 20 #define GCTAINSTDIR_HPP
21 21  
22 22 /* __ Includes ___________________________________________________________ */
23   -#include <iostream>
  23 +#include <string>
24 24 #include "GInstDir.hpp"
25 25 #include "GSkyDir.hpp"
26 26  
... ... @@ -35,9 +35,6 @@
35 35 ***************************************************************************/
36 36 class GCTAInstDir : public GInstDir {
37 37  
38   - // I/O friends
39   - friend std::ostream& operator<< (std::ostream& os, const GCTAInstDir& dir);
40   -
41 38 public:
42 39 // Constructors and destructors
43 40 GCTAInstDir(void);
... ... @@ -48,36 +45,39 @@ public:
48 45 // Operators
49 46 GCTAInstDir& operator= (const GCTAInstDir& dir);
50 47  
51   - // Methods
52   - void clear(void);
53   - void skydir(const GSkyDir& dir) { m_dir=dir; }
54   - void radec(const double& ra, const double& dec) { m_dir.radec(ra,dec); }
55   - void radec_deg(const double& ra, const double& dec) { m_dir.radec_deg(ra,dec); }
56   - void lb(const double& l, const double& b) { m_dir.lb(l,b); }
57   - void lb_deg(const double& l, const double& b) { m_dir.lb_deg(l,b); }
58   - GSkyDir skydir(void) const { return m_dir; }
59   - double l(void) const { return m_dir.l(); }
60   - double l_deg(void) const { return m_dir.l_deg(); }
61   - double b(void) const { return m_dir.b(); }
62   - double b_deg(void) const { return m_dir.b_deg(); }
63   - double ra(void) const { return m_dir.ra(); }
64   - double ra_deg(void) const { return m_dir.ra_deg(); }
65   - double dec(void) const { return m_dir.dec(); }
66   - double dec_deg(void) const { return m_dir.dec_deg(); }
67   - double dist(const GSkyDir& dir) const { return m_dir.dist(dir); }
68   - double dist_deg(const GSkyDir& dir) const { return m_dir.dist_deg(dir); }
69   - double dist(const GCTAInstDir& dir) const;
70   - double dist_deg(const GCTAInstDir& dir) const;
  48 + // Implemented pure virtual base class methods
  49 + void clear(void);
  50 + GCTAInstDir* clone(void) const;
  51 + std::string print(void) const;
  52 +
  53 + // Other methods
  54 + void skydir(const GSkyDir& dir) { m_dir=dir; }
  55 + void radec(const double& ra, const double& dec) { m_dir.radec(ra,dec); }
  56 + void radec_deg(const double& ra, const double& dec) { m_dir.radec_deg(ra,dec); }
  57 + void lb(const double& l, const double& b) { m_dir.lb(l,b); }
  58 + void lb_deg(const double& l, const double& b) { m_dir.lb_deg(l,b); }
  59 + GSkyDir skydir(void) const { return m_dir; }
  60 + double l(void) const { return m_dir.l(); }
  61 + double l_deg(void) const { return m_dir.l_deg(); }
  62 + double b(void) const { return m_dir.b(); }
  63 + double b_deg(void) const { return m_dir.b_deg(); }
  64 + double ra(void) const { return m_dir.ra(); }
  65 + double ra_deg(void) const { return m_dir.ra_deg(); }
  66 + double dec(void) const { return m_dir.dec(); }
  67 + double dec_deg(void) const { return m_dir.dec_deg(); }
  68 + double dist(const GSkyDir& dir) const { return m_dir.dist(dir); }
  69 + double dist_deg(const GSkyDir& dir) const { return m_dir.dist_deg(dir); }
  70 + double dist(const GCTAInstDir& dir) const;
  71 + double dist_deg(const GCTAInstDir& dir) const;
71 72  
72 73 protected:
73 74 // Protected methods
74   - void init_members(void);
75   - void copy_members(const GCTAInstDir& dir);
76   - void free_members(void);
77   - virtual GCTAInstDir* clone(void) const;
  75 + void init_members(void);
  76 + void copy_members(const GCTAInstDir& dir);
  77 + void free_members(void);
78 78  
79 79 // Data members
80   - GSkyDir m_dir; //!< Incident direction of event
  80 + GSkyDir m_dir; //!< Observed incident direction of event
81 81 };
82 82  
83 83 #endif /* GCTAINSTDIR_HPP */
... ...
inst/cta/pyext/ctalib.i
  View file @5fa82e7
... ... @@ -16,9 +16,11 @@
16 16 %module ctalib
17 17  
18 18 /* __ CTA ________________________________________________________________ */
  19 +#include "GCTAObservation.i"
19 20 #include "GCTAPointing.i"
20 21 #include "GCTAResponse.i"
21   -#include "GCTAObservation.i"
  22 +#include "GCTAInstDir.i"
  23 +#include "GCTARoi.i"
22 24  
23 25  
24 26  
... ...
inst/cta/test/test_CTA.cpp
  View file @5fa82e7
... ... @@ -52,7 +52,7 @@ GModels crab_plaw(void)
52 52 dir.radec_deg(117.0, -33.0); // Adapt to source position in file
53 53 point_source = GModelSpatialPtsrc(dir);
54 54 power_law = GModelSpectralPlaw(1.0e-7, -2.1);
55   - power_law.par(0)->min(1.0e-12);
  55 + power_law(0).min(1.0e-12);
56 56 crab = GModel(point_source, power_law);
57 57 crab.name("Crab");
58 58 models.append(crab);
... ...
inst/lat/Makefile.am
  View file @5fa82e7
... ... @@ -19,7 +19,8 @@ sources = src/GLATObservation.cpp \
19 19 src/GLATResponseTable.cpp \
20 20 src/GLATAeff.cpp \
21 21 src/GLATPsf.cpp \
22   - src/GLATEdisp.cpp
  22 + src/GLATEdisp.cpp \
  23 + src/GLATException.cpp
23 24  
24 25 # Define headers to be installed
25 26 pkginclude_HEADERS = include/GLATObservation.hpp \
... ... @@ -32,7 +33,8 @@ pkginclude_HEADERS = include/GLATObservation.hpp \
32 33 include/GLATRoi.hpp \
33 34 include/GLATResponse.hpp \
34 35 include/GLATResponseTable.hpp \
35   - include/GLATLib.hpp
  36 + include/GLATLib.hpp \
  37 + include/GLATException.hpp
36 38  
37 39 # Build libtool library
38 40 noinst_LTLIBRARIES = liblat.la
... ...
inst/lat/include/GLATInstDir.hpp
  View file @5fa82e7
... ... @@ -20,7 +20,7 @@
20 20 #define GLATINSTDIR_HPP
21 21  
22 22 /* __ Includes ___________________________________________________________ */
23   -#include <iostream>
  23 +#include <string>
24 24 #include "GInstDir.hpp"
25 25 #include "GSkyDir.hpp"
26 26  
... ... @@ -35,46 +35,49 @@
35 35 ***************************************************************************/
36 36 class GLATInstDir : public GInstDir {
37 37  
38   - // I/O friends
39   - friend std::ostream& operator<< (std::ostream& os, const GLATInstDir& dir);
40   -
41 38 public:
42 39 // Constructors and destructors
43 40 GLATInstDir(void);
  41 + explicit GLATInstDir(const GSkyDir& dir);
44 42 GLATInstDir(const GLATInstDir& dir);
45 43 virtual ~GLATInstDir(void);
46 44  
47 45 // Operators
48 46 GLATInstDir& operator= (const GLATInstDir& dir);
49 47  
  48 + // Implemented pure virtual base class methods
  49 + void clear(void);
  50 + GLATInstDir* clone(void) const;
  51 + std::string print(void) const;
  52 +
50 53 // Methods
51   - void clear(void);
52   - void radec(const double& ra, const double& dec) { m_dir.radec(ra,dec); }
53   - void radec_deg(const double& ra, const double& dec) { m_dir.radec_deg(ra,dec); }
54   - void lb(const double& l, const double& b) { m_dir.lb(l,b); }
55   - void lb_deg(const double& l, const double& b) { m_dir.lb_deg(l,b); }
56   - double l(void) const { return m_dir.l(); }
57   - double l_deg(void) const { return m_dir.l_deg(); }
58   - double b(void) const { return m_dir.b(); }
59   - double b_deg(void) const { return m_dir.b_deg(); }
60   - double ra(void) const { return m_dir.ra(); }
61   - double ra_deg(void) const { return m_dir.ra_deg(); }
62   - double dec(void) const { return m_dir.dec(); }
63   - double dec_deg(void) const { return m_dir.dec_deg(); }
64   - double dist(GSkyDir& dir) const { return m_dir.dist(dir); }
65   - double dist_deg(GSkyDir& dir) const { return m_dir.dist_deg(dir); }
66   - double dist(GLATInstDir& dir) const;
67   - double dist_deg(GLATInstDir& dir) const;
  54 + void skydir(const GSkyDir& dir) { m_dir=dir; }
  55 + void radec(const double& ra, const double& dec) { m_dir.radec(ra,dec); }
  56 + void radec_deg(const double& ra, const double& dec) { m_dir.radec_deg(ra,dec); }
  57 + void lb(const double& l, const double& b) { m_dir.lb(l,b); }
  58 + void lb_deg(const double& l, const double& b) { m_dir.lb_deg(l,b); }
  59 + GSkyDir skydir(void) const { return m_dir; }
  60 + double l(void) const { return m_dir.l(); }
  61 + double l_deg(void) const { return m_dir.l_deg(); }
  62 + double b(void) const { return m_dir.b(); }
  63 + double b_deg(void) const { return m_dir.b_deg(); }
  64 + double ra(void) const { return m_dir.ra(); }
  65 + double ra_deg(void) const { return m_dir.ra_deg(); }
  66 + double dec(void) const { return m_dir.dec(); }
  67 + double dec_deg(void) const { return m_dir.dec_deg(); }
  68 + double dist(GSkyDir& dir) const { return m_dir.dist(dir); }
  69 + double dist_deg(GSkyDir& dir) const { return m_dir.dist_deg(dir); }
  70 + double dist(GLATInstDir& dir) const;
  71 + double dist_deg(GLATInstDir& dir) const;
68 72  
69 73 protected:
70 74 // Protected methods
71   - void init_members(void);
72   - void copy_members(const GLATInstDir& dir);
73   - void free_members(void);
74   - virtual GLATInstDir* clone(void) const;
  75 + void init_members(void);
  76 + void copy_members(const GLATInstDir& dir);
  77 + void free_members(void);
75 78  
76 79 // Data members
77   - GSkyDir m_dir; //!< Incident direction of event
  80 + GSkyDir m_dir; //!< Observed incident direction of event
78 81 };
79 82  
80 83 #endif /* GLATINSTDIR_HPP */
... ...
inst/lat/pyext/latlib.i
  View file @5fa82e7
... ... @@ -17,7 +17,11 @@
17 17  
18 18 /* __ LAT ________________________________________________________________ */
19 19 #include "GLATObservation.i"
20   -//#include "GLATResponse.i"
  20 +#include "GLATEventCube.i"
  21 +#include "GLATEventBin.i"
  22 +#include "GLATResponse.i"
  23 +#include "GLATInstDir.i"
  24 +#include "GLATRoi.i"
21 25  
22 26  
23 27  
... ...
inst/lat/src/GLATInstDir.cpp
  View file @5fa82e7
... ... @@ -20,8 +20,8 @@
20 20 #ifdef HAVE_CONFIG_H
21 21 #include <config.h>
22 22 #endif
23   -#include "GTools.hpp"
24 23 #include "GLATInstDir.hpp"
  24 +#include "GTools.hpp"
25 25  
26 26 /* __ Method name definitions ____________________________________________ */
27 27  
... ... @@ -40,9 +40,9 @@
40 40 ==========================================================================*/
41 41  
42 42 /***********************************************************************//**
43   - * @brief Constructor
  43 + * @brief Void constructor
44 44 ***************************************************************************/
45   -GLATInstDir::GLATInstDir(void)
  45 +GLATInstDir::GLATInstDir(void) : GInstDir()
46 46 {
47 47 // Initialise class members
48 48 init_members();
... ... @@ -53,12 +53,31 @@ GLATInstDir::GLATInstDir(void)
53 53  
54 54  
55 55 /***********************************************************************//**
  56 + * @brief GSkyDir constructor
  57 + *
  58 + * @param[in] dir Sky direction.
  59 + *
  60 + * Construct LAT instrument direction from sky direction.
  61 + ***************************************************************************/
  62 +GLATInstDir::GLATInstDir(const GSkyDir& dir) : GInstDir()
  63 +{
  64 + // Initialise class members
  65 + init_members();
  66 +
  67 + // Assign sky direction
  68 + m_dir = dir;
  69 +
  70 + // Return
  71 + return;
  72 +}
  73 +
  74 +
  75 +/***********************************************************************//**
56 76 * @brief Copy constructor
57 77 *
58   - * @param[in] dir Instrument direction from which class should be
59   - * instantiated.
  78 + * @param[in] dir Instrument direction.
60 79 ***************************************************************************/
61   -GLATInstDir::GLATInstDir(const GLATInstDir& dir)
  80 +GLATInstDir::GLATInstDir(const GLATInstDir& dir) : GInstDir(dir)
62 81 {
63 82 // Initialise class members
64 83 init_members();
... ... @@ -93,13 +112,16 @@ GLATInstDir::~GLATInstDir(void)
93 112 /***********************************************************************//**
94 113 * @brief Assignment operator
95 114 *
96   - * @param[in] dir Instrument direction to be assigned.
  115 + * @param[in] dir Instrument direction.
97 116 ***************************************************************************/
98 117 GLATInstDir& GLATInstDir::operator= (const GLATInstDir& dir)
99 118 {
100 119 // Execute only if object is not identical
101 120 if (this != &dir) {
102 121  
  122 + // Copy base class members
  123 + this->GInstDir::operator=(dir);
  124 +
103 125 // Free members
104 126 free_members();
105 127  
... ... @@ -123,14 +145,16 @@ GLATInstDir&amp; GLATInstDir::operator= (const GLATInstDir&amp; dir)
123 145 ==========================================================================*/
124 146  
125 147 /***********************************************************************//**
126   - * @brief Clear sky direction
  148 + * @brief Clear instance
127 149 ***************************************************************************/
128 150 void GLATInstDir::clear(void)
129 151 {
130 152 // Free members
131 153 free_members();
  154 + this->GInstDir::free_members();
132 155  
133 156 // Initialise private members
  157 + this->GInstDir::init_members();
134 158 init_members();
135 159  
136 160 // Return
... ... @@ -139,9 +163,18 @@ void GLATInstDir::clear(void)
139 163  
140 164  
141 165 /***********************************************************************//**
  166 + * @brief Clone instance
  167 + ***************************************************************************/
  168 +GLATInstDir* GLATInstDir::clone(void) const
  169 +{
  170 + return new GLATInstDir(*this);
  171 +}
  172 +
  173 +
  174 +/***********************************************************************//**
142 175 * @brief Compute angular distance between instrument directions in radians
143 176 *
144   - * @param[in] dir Instrument direction to which distance is to be computed.
  177 + * @param[in] dir Instrument direction.
145 178 ***************************************************************************/
146 179 double GLATInstDir::dist(GLATInstDir& dir) const
147 180 {
... ... @@ -162,7 +195,7 @@ double GLATInstDir::dist(GLATInstDir&amp; dir) const
162 195 /***********************************************************************//**
163 196 * @brief Compute angular distance between instrument directions in degrees
164 197 *
165   - * @param[in] dir Instrument direction to which distance is to be computed.
  198 + * @param[in] dir Instrument direction.
166 199 ***************************************************************************/
167 200 double GLATInstDir::dist_deg(GLATInstDir& dir) const
168 201 {
... ... @@ -171,6 +204,22 @@ double GLATInstDir::dist_deg(GLATInstDir&amp; dir) const
171 204 }
172 205  
173 206  
  207 +/***********************************************************************//**
  208 + * @brief Print instrument direction information
  209 + ***************************************************************************/
  210 +std::string GLATInstDir::print(void) const
  211 +{
  212 + // Initialise result string
  213 + std::string result;
  214 +
  215 + // Append instrument direction
  216 + result.append("RA="+str(ra_deg())+", DEC="+str(dec_deg()));
  217 +
  218 + // Return result
  219 + return result;
  220 +}
  221 +
  222 +
174 223 /*==========================================================================
175 224 = =
176 225 = Private methods =
... ... @@ -193,7 +242,7 @@ void GLATInstDir::init_members(void)
193 242 /***********************************************************************//**
194 243 * @brief Copy class members
195 244 *
196   - * @param[in] dir Instrument direction from which members should be copied
  245 + * @param[in] dir Instrument direction.
197 246 ***************************************************************************/
198 247 void GLATInstDir::copy_members(const GLATInstDir& dir)
199 248 {
... ... @@ -215,32 +264,8 @@ void GLATInstDir::free_members(void)
215 264 }
216 265  
217 266  
218   -/***********************************************************************//**
219   - * @brief Clone instance
220   - ***************************************************************************/
221   -GLATInstDir* GLATInstDir::clone(void) const
222   -{
223   - return new GLATInstDir(*this);
224   -}
225   -
226   -
227 267 /*==========================================================================
228 268 = =
229 269 = Friends =
230 270 = =
231 271 ==========================================================================*/
232   -
233   -/***********************************************************************//**
234   - * @brief Output operator
235   - *
236   - * @param[in] os Output stream
237   - * @param[in] column Instrument direction to put in output stream
238   - ***************************************************************************/
239   -std::ostream& operator<< (std::ostream& os, const GLATInstDir& dir)
240   -{
241   - // Output sky direction
242   - os << dir.m_dir;
243   -
244   - // Return output stream
245   - return os;
246   -}
... ...
inst/mwl/test/test_python.py
  View file @5fa82e7
... ... @@ -87,6 +87,7 @@ def fit_spectrum(filename, xmlname):
87 87 x = np.power(10., np.arange(-1., 8., 0.1))
88 88 y = []
89 89 d = GSkyDir()
  90 + d.radec_deg(83.6331, 22.0145)
90 91 e = GEnergy()
91 92 t = GTime()
92 93 for energy in x:
... ... @@ -121,7 +122,8 @@ if __name__ == &#39;__main__&#39;:
121 122 extract_data("data/crab_mwl.fits")
122 123  
123 124 # Fit spectrum
124   - fit_spectrum("data/crab_mwl.fits", "data/crab_mwl2.xml")
  125 +# fit_spectrum("data/crab_mwl.fits", "data/crab_mwl2.xml")
  126 + fit_spectrum("data/crab_mwl.fits", "data/crab_mwl3.xml")
125 127  
126 128 # Show spectrum
127 129 plt.show()
... ...
pyext/GOptimizerLM.i
  View file @5fa82e7
... ... @@ -31,6 +31,7 @@ public:
31 31  
32 32 // Constructors and destructors
33 33 GOptimizerLM(void);
  34 + GOptimizerLM(GLog& log);
34 35 GOptimizerLM(const GOptimizerLM& opt);
35 36 virtual ~GOptimizerLM(void);
36 37  
... ...
pyext/gammalib.i
  View file @5fa82e7
... ... @@ -64,7 +64,7 @@
64 64  
65 65 /* __ Application module _________________________________________________ */
66 66 #include "GApplication.i"
67   -//#include "GLog.i"
  67 +#include "GLog.i"
68 68 #include "GPars.i"
69 69 #include "GPar.i"
70 70  
... ... @@ -113,6 +113,7 @@
113 113 #include "GModelSpatialPtsrc.i"
114 114 #include "GModelSpectral.i"
115 115 #include "GModelSpectralPlaw.i"
  116 +#include "GModelSpectralExpPlaw.i"
116 117 #include "GModelTemporal.i"
117 118 #include "GModelTemporalConst.i"
118 119  
... ...
src/obs/GObservations_optimizer.cpp
  View file @5fa82e7
... ... @@ -34,8 +34,9 @@
34 34 /* __ Coding definitions _________________________________________________ */
35 35  
36 36 /* __ Debug definitions __________________________________________________ */
37   -#define G_EVAL_TIMING 0 //!< Perform optimizer timing (0=no, 1=yes)
38   -#define G_EVAL_DEBUG 0 //!< Perform optimizer debugging (0=no, 1=yes)
  37 +#define G_EVAL_TIMING 0 //!< Perform optimizer timing (0=no, 1=yes)
  38 +#define G_EVAL_DEBUG 0 //!< Perform optimizer debugging (0=no, 1=yes)
  39 +#define G_OPT_DEBUG 0 //!< Perform optimizer debugging (0=no, 1=yes)
39 40  
40 41 /* __ Prototypes _________________________________________________________ */
41 42  
... ... @@ -47,7 +48,7 @@
47 48 ==========================================================================*/
48 49  
49 50 /***********************************************************************//**
50   - * @brief Constructor
  51 + * @brief Void constructor
51 52 ***************************************************************************/
52 53 GObservations::optimizer::optimizer(void) : GOptimizerFunction()
53 54 {
... ... @@ -99,7 +100,7 @@ GObservations::optimizer::optimizer(const optimizer&amp; fct) : GOptimizerFunction(f
99 100 /***********************************************************************//**
100 101 * @brief Destructor
101 102 ***************************************************************************/
102   -GObservations::optimizer::~optimizer()
  103 +GObservations::optimizer::~optimizer(void)
103 104 {
104 105 // Free members
105 106 free_members();
... ... @@ -118,7 +119,7 @@ GObservations::optimizer::~optimizer()
118 119 /***********************************************************************//**
119 120 * @brief Assignment operator
120 121 *
121   - * @param[in] obs Instance to be assigned
  122 + * @param[in] obs Instance to be assigned.
122 123 ***************************************************************************/
123 124 GObservations::optimizer& GObservations::optimizer::operator= (const optimizer& fct)
124 125 {
... ... @@ -209,13 +210,17 @@ void GObservations::optimizer::eval(const GOptimizerPars&amp; pars)
209 210 m_npred += m_this->m_obs[i]->npred((GModels&)pars, m_wrk_grad);
210 211 *m_gradient += *m_wrk_grad;
211 212 #if G_EVAL_DEBUG
212   - std::cout << "Npred=" << m_npred << " Grad="
213   - << *m_gradient << std::endl;
  213 + std::cout << "Unbinned Poisson";
  214 + std::cout << " Npred=" << m_npred;
  215 + std::cout << " Grad="<< *m_gradient << std::endl;
214 216 #endif
215 217  
216 218 // Update the log-likelihood
217 219 poisson_unbinned(*(m_this->m_obs[i]), pars);
218 220  
  221 + // Add the Npred value to the log-likelihood
  222 + m_value += m_npred;
  223 +
219 224 } // endif: Poisson statistics
220 225  
221 226 // ... otherwise throw an exception
... ... @@ -229,12 +234,20 @@ void GObservations::optimizer::eval(const GOptimizerPars&amp; pars)
229 234 else {
230 235  
231 236 // Poisson statistics
232   - if (toupper(statistics) == "POISSON")
  237 + if (toupper(statistics) == "POISSON") {
  238 + #if G_EVAL_DEBUG
  239 + std::cout << "Binned Poisson" << std::endl;
  240 + #endif
233 241 poisson_binned(*(m_this->m_obs[i]), pars);
  242 + }
234 243  
235 244 // ... or Gaussian statistics
236   - else if (toupper(statistics) == "GAUSSIAN")
  245 + else if (toupper(statistics) == "GAUSSIAN") {
  246 + #if G_EVAL_DEBUG
  247 + std::cout << "Binned Gaussian" << std::endl;
  248 + #endif
237 249 gaussian_binned(*(m_this->m_obs[i]), pars);
  250 + }
238 251  
239 252 // ... or unsupported
240 253 else
... ... @@ -245,9 +258,6 @@ void GObservations::optimizer::eval(const GOptimizerPars&amp; pars)
245 258  
246 259 } // endfor: looped over observations
247 260  
248   - // Add the Npred value to the log-likelihood
249   - m_value += m_npred;
250   -
251 261 // Release stack
252 262 m_covar->stack_destroy();
253 263  
... ... @@ -311,9 +321,7 @@ void GObservations::optimizer::poisson_unbinned(const GObservation&amp; obs,
311 321 GEvent* event = obs.events()->pointer(i);
312 322  
313 323 // Get model and derivative
314   - double model = obs.model((GModels&)pars,
315   - event->dir(), event->energy(), event->time(),
316   - m_wrk_grad);
  324 + double model = obs.model((GModels&)pars, *event, m_wrk_grad);
317 325  
318 326 // Skip bin if model is too small (avoids -Inf or NaN gradients)
319 327 if (model <= m_minmod)
... ... @@ -410,6 +418,17 @@ void GObservations::optimizer::poisson_binned(const GObservation&amp; obs,
410 418 clock_t t_start = clock();
411 419 #endif
412 420  
  421 + // Initialise statistics
  422 + #if G_OPT_DEBUG
  423 + int n_bins = 0;
  424 + int n_used = 0;
  425 + int n_small_model = 0;
  426 + int n_zero_data = 0;
  427 + double sum_data = 0.0;
  428 + double sum_model = 0.0;
  429 + double init_value = m_value;
  430 + #endif
  431 +
413 432 // Get number of parameters
414 433 int npars = pars.npars();
415 434  
... ... @@ -420,6 +439,11 @@ void GObservations::optimizer::poisson_binned(const GObservation&amp; obs,
420 439 // Iterate over all bins
421 440 for (int i = 0; i < obs.events()->size(); ++i) {
422 441  
  442 + // Update number of bins
  443 + #if G_OPT_DEBUG
  444 + n_bins++;
  445 + #endif
  446 +
423 447 // Get pointer to bin
424 448 GEventBin* bin = (GEventBin*)obs.events()->pointer(i);
425 449  
... ... @@ -427,16 +451,25 @@ void GObservations::optimizer::poisson_binned(const GObservation&amp; obs,
427 451 double data = bin->counts();
428 452  
429 453 // Get model and derivative
430   - double model = obs.model((GModels&)pars,
431   - bin->dir(), bin->energy(), bin->time(),
432   - m_wrk_grad);
  454 + double model = obs.model((GModels&)pars, *bin, m_wrk_grad);
433 455  
434 456 // Multiply model by bin size
435 457 model *= bin->size();
436 458  
437 459 // Skip bin if model is too small (avoids -Inf or NaN gradients)
438   - if (model <= m_minmod)
  460 + if (model <= m_minmod) {
  461 + #if G_OPT_DEBUG
  462 + n_small_model++;
  463 + #endif
439 464 continue;
  465 + }
  466 +
  467 + // Update statistics
  468 + #if G_OPT_DEBUG
  469 + n_used++;
  470 + sum_data += data;
  471 + sum_model += model;
  472 + #endif
440 473  
441 474 // Update Npred
442 475 m_npred += model;
... ... @@ -498,6 +531,11 @@ void GObservations::optimizer::poisson_binned(const GObservation&amp; obs,
498 531 // ... handle now data=0
499 532 else {
500 533  
  534 + // Update statistics
  535 + #if G_OPT_DEBUG
  536 + n_zero_data++;
  537 + #endif
  538 +
501 539 // Update Poissonian statistics (excluding factorial term for
502 540 // faster computation)
503 541 m_value += model;
... ... @@ -519,6 +557,18 @@ void GObservations::optimizer::poisson_binned(const GObservation&amp; obs,
519 557 if (values != NULL) delete [] values;
520 558 if (inx != NULL) delete [] inx;
521 559  
  560 + // Dump statistics
  561 + #if G_OPT_DEBUG
  562 + std::cout << "Number of bins: " << n_bins << std::endl;
  563 + std::cout << "Number of bins used for computation: " << n_used << std::endl;
  564 + std::cout << "Number of bins excluded due to small model: " << n_small_model << std::endl;
  565 + std::cout << "Number of bins with zero data: " << n_zero_data << std::endl;
  566 + std::cout << "Sum of data: " << sum_data << std::endl;
  567 + std::cout << "Sum of model: " << sum_model << std::endl;
  568 + std::cout << "Initial statistics: " << init_value << std::endl;
  569 + std::cout << "Statistics: " << m_value-init_value << std::endl;
  570 + #endif
  571 +
522 572 // Optionally dump gradient and covariance matrix
523 573 #if G_EVAL_DEBUG
524 574 std::cout << *m_gradient << std::endl;
... ... @@ -589,9 +639,7 @@ void GObservations::optimizer::gaussian_binned(const GObservation&amp; obs,
589 639 continue;
590 640  
591 641 // Get model and derivative
592   - double model = obs.model((GModels&)pars,
593   - bin->dir(), bin->energy(), bin->time(),
594   - m_wrk_grad);
  642 + double model = obs.model((GModels&)pars, *bin, m_wrk_grad);
595 643  
596 644 // Multiply model by bin size
597 645 model *= bin->size();
... ...