diff --git a/include/GEvent.hpp b/include/GEvent.hpp index ccc1373..0972319 100644 --- a/include/GEvent.hpp +++ b/include/GEvent.hpp @@ -20,6 +20,9 @@ #define GEVENT_HPP /* __ Includes ___________________________________________________________ */ +#include +#include +#include "GLog.hpp" #include "GInstDir.hpp" #include "GEnergy.hpp" #include "GTime.hpp" @@ -60,6 +63,10 @@ ***************************************************************************/ class GEvent { + // I/O friends + friend std::ostream& operator<< (std::ostream& os, const GEvent& event); + friend GLog& operator<< (GLog& log, const GEvent& event); + public: // Constructors and destructors GEvent(void); @@ -69,19 +76,18 @@ public: // Operators virtual GEvent& operator= (const GEvent& event); - // Event access methods + // Pure virtual methods + virtual void clear(void) = 0; + virtual GEvent* clone(void) const = 0; + virtual double size(void) const = 0; virtual const GInstDir& dir(void) const = 0; virtual const GEnergy& energy(void) const = 0; virtual const GTime& time(void) const = 0; virtual double counts(void) const = 0; virtual double error(void) const = 0; - - // Other methods - virtual void clear(void) = 0; - virtual GEvent* clone(void) const = 0; - virtual double size(void) const = 0; - virtual bool isatom(void) const = 0; - virtual bool isbin(void) const = 0; + virtual bool isatom(void) const = 0; + virtual bool isbin(void) const = 0; + virtual std::string print(void) const = 0; protected: // Protected methods diff --git a/include/GEventAtom.hpp b/include/GEventAtom.hpp index 0c49505..816a537 100644 --- a/include/GEventAtom.hpp +++ b/include/GEventAtom.hpp @@ -20,6 +20,7 @@ #define GEVENTATOM_HPP /* __ Includes ___________________________________________________________ */ +#include #include "GEvent.hpp" #include "GInstDir.hpp" #include "GEnergy.hpp" @@ -64,19 +65,20 @@ public: // Operators virtual GEventAtom& operator= (const GEventAtom& atom); - // Event access methods + // Pure virtual methods + virtual void clear(void) = 0; + virtual GEvent* clone(void) const = 0; + virtual double size(void) const { return 1.0; } virtual const GInstDir& dir(void) const = 0; virtual const GEnergy& energy(void) const = 0; virtual const GTime& time(void) const = 0; - double counts(void) const { return 1.0; } - double error(void) const { return 0.0; } + virtual double counts(void) const { return 1.0; } + virtual double error(void) const { return 0.0; } + virtual std::string print(void) const = 0; // Other methods - virtual void clear(void) = 0; - virtual double size(void) const { return 1.0; } - virtual GEventAtom* clone(void) const = 0; - virtual bool isatom(void) const { return true; } - virtual bool isbin(void) const { return false; } + bool isatom(void) const { return true; } + bool isbin(void) const { return false; } protected: // Protected methods diff --git a/include/GEventBin.hpp b/include/GEventBin.hpp index 3b84bcc..800d8fd 100644 --- a/include/GEventBin.hpp +++ b/include/GEventBin.hpp @@ -20,6 +20,7 @@ #define GEVENTBIN_HPP /* __ Includes ___________________________________________________________ */ +#include #include "GEvent.hpp" #include "GInstDir.hpp" #include "GEnergy.hpp" @@ -66,19 +67,20 @@ public: // Operators virtual GEventBin& operator= (const GEventBin& bin); - // Event access methods + // Pure virtual methods + virtual void clear(void) = 0; + virtual GEvent* clone(void) const = 0; + virtual double size(void) const = 0; virtual const GInstDir& dir(void) const = 0; virtual const GEnergy& energy(void) const = 0; virtual const GTime& time(void) const = 0; virtual double counts(void) const = 0; - virtual double error(void) const; + virtual double error(void) const = 0; + virtual std::string print(void) const = 0; // Other methods - virtual void clear(void) = 0; - virtual double size(void) const = 0; - virtual GEventBin* clone(void) const = 0; - bool isatom(void) const { return false; } - bool isbin(void) const { return true; } + bool isatom(void) const { return false; } + bool isbin(void) const { return true; } protected: // Protected methods diff --git a/include/GammaLib.hpp b/include/GammaLib.hpp index 0c5d36b..7e0e00d 100644 --- a/include/GammaLib.hpp +++ b/include/GammaLib.hpp @@ -124,9 +124,14 @@ #include "GModelPar.hpp" #include "GModelSpatial.hpp" #include "GModelSpatialPtsrc.hpp" +#include "GModelSpatialConst.hpp" +#include "GModelSpatialCube.hpp" #include "GModelSpectral.hpp" #include "GModelSpectralPlaw.hpp" +#include "GModelSpectralPlaw2.hpp" #include "GModelSpectralExpPlaw.hpp" +#include "GModelSpectralFunc.hpp" +#include "GModelSpectralConst.hpp" #include "GModelTemporal.hpp" #include "GModelTemporalConst.hpp" diff --git a/include/Makefile.am b/include/Makefile.am index f73e872..d9eea39 100644 --- a/include/Makefile.am +++ b/include/Makefile.am @@ -29,9 +29,14 @@ pkginclude_HEADERS = GammaLib.hpp \ GModelPar.hpp \ GModelSpatial.hpp \ GModelSpatialPtsrc.hpp \ + GModelSpatialConst.hpp \ + GModelSpatialCube.hpp \ GModelSpectral.hpp \ GModelSpectralPlaw.hpp \ + GModelSpectralPlaw2.hpp \ GModelSpectralExpPlaw.hpp \ + GModelSpectralFunc.hpp \ + GModelSpectralConst.hpp \ GModelTemporal.hpp \ GModelTemporalConst.hpp \ GApplication.hpp \ diff --git a/inst/cta/include/GCTAInstDir.hpp b/inst/cta/include/GCTAInstDir.hpp index b5f6860..2b6193b 100644 --- a/inst/cta/include/GCTAInstDir.hpp +++ b/inst/cta/include/GCTAInstDir.hpp @@ -20,7 +20,7 @@ #define GCTAINSTDIR_HPP /* __ Includes ___________________________________________________________ */ -#include +#include #include "GInstDir.hpp" #include "GSkyDir.hpp" @@ -35,9 +35,6 @@ ***************************************************************************/ class GCTAInstDir : public GInstDir { - // I/O friends - friend std::ostream& operator<< (std::ostream& os, const GCTAInstDir& dir); - public: // Constructors and destructors GCTAInstDir(void); @@ -48,36 +45,39 @@ public: // Operators GCTAInstDir& operator= (const GCTAInstDir& dir); - // Methods - void clear(void); - void skydir(const GSkyDir& dir) { m_dir=dir; } - void radec(const double& ra, const double& dec) { m_dir.radec(ra,dec); } - void radec_deg(const double& ra, const double& dec) { m_dir.radec_deg(ra,dec); } - void lb(const double& l, const double& b) { m_dir.lb(l,b); } - void lb_deg(const double& l, const double& b) { m_dir.lb_deg(l,b); } - GSkyDir skydir(void) const { return m_dir; } - double l(void) const { return m_dir.l(); } - double l_deg(void) const { return m_dir.l_deg(); } - double b(void) const { return m_dir.b(); } - double b_deg(void) const { return m_dir.b_deg(); } - double ra(void) const { return m_dir.ra(); } - double ra_deg(void) const { return m_dir.ra_deg(); } - double dec(void) const { return m_dir.dec(); } - double dec_deg(void) const { return m_dir.dec_deg(); } - double dist(const GSkyDir& dir) const { return m_dir.dist(dir); } - double dist_deg(const GSkyDir& dir) const { return m_dir.dist_deg(dir); } - double dist(const GCTAInstDir& dir) const; - double dist_deg(const GCTAInstDir& dir) const; + // Implemented pure virtual base class methods + void clear(void); + GCTAInstDir* clone(void) const; + std::string print(void) const; + + // Other methods + void skydir(const GSkyDir& dir) { m_dir=dir; } + void radec(const double& ra, const double& dec) { m_dir.radec(ra,dec); } + void radec_deg(const double& ra, const double& dec) { m_dir.radec_deg(ra,dec); } + void lb(const double& l, const double& b) { m_dir.lb(l,b); } + void lb_deg(const double& l, const double& b) { m_dir.lb_deg(l,b); } + GSkyDir skydir(void) const { return m_dir; } + double l(void) const { return m_dir.l(); } + double l_deg(void) const { return m_dir.l_deg(); } + double b(void) const { return m_dir.b(); } + double b_deg(void) const { return m_dir.b_deg(); } + double ra(void) const { return m_dir.ra(); } + double ra_deg(void) const { return m_dir.ra_deg(); } + double dec(void) const { return m_dir.dec(); } + double dec_deg(void) const { return m_dir.dec_deg(); } + double dist(const GSkyDir& dir) const { return m_dir.dist(dir); } + double dist_deg(const GSkyDir& dir) const { return m_dir.dist_deg(dir); } + double dist(const GCTAInstDir& dir) const; + double dist_deg(const GCTAInstDir& dir) const; protected: // Protected methods - void init_members(void); - void copy_members(const GCTAInstDir& dir); - void free_members(void); - virtual GCTAInstDir* clone(void) const; + void init_members(void); + void copy_members(const GCTAInstDir& dir); + void free_members(void); // Data members - GSkyDir m_dir; //!< Incident direction of event + GSkyDir m_dir; //!< Observed incident direction of event }; #endif /* GCTAINSTDIR_HPP */ diff --git a/inst/cta/pyext/ctalib.i b/inst/cta/pyext/ctalib.i index 7060f93..3cc1db3 100755 --- a/inst/cta/pyext/ctalib.i +++ b/inst/cta/pyext/ctalib.i @@ -16,9 +16,11 @@ %module ctalib /* __ CTA ________________________________________________________________ */ +#include "GCTAObservation.i" #include "GCTAPointing.i" #include "GCTAResponse.i" -#include "GCTAObservation.i" +#include "GCTAInstDir.i" +#include "GCTARoi.i" diff --git a/inst/cta/test/test_CTA.cpp b/inst/cta/test/test_CTA.cpp index 28b2f58..12600ad 100644 --- a/inst/cta/test/test_CTA.cpp +++ b/inst/cta/test/test_CTA.cpp @@ -52,7 +52,7 @@ GModels crab_plaw(void) dir.radec_deg(117.0, -33.0); // Adapt to source position in file point_source = GModelSpatialPtsrc(dir); power_law = GModelSpectralPlaw(1.0e-7, -2.1); - power_law.par(0)->min(1.0e-12); + power_law(0).min(1.0e-12); crab = GModel(point_source, power_law); crab.name("Crab"); models.append(crab); diff --git a/inst/lat/Makefile.am b/inst/lat/Makefile.am index 00805d9..6bf3c2c 100644 --- a/inst/lat/Makefile.am +++ b/inst/lat/Makefile.am @@ -19,7 +19,8 @@ sources = src/GLATObservation.cpp \ src/GLATResponseTable.cpp \ src/GLATAeff.cpp \ src/GLATPsf.cpp \ - src/GLATEdisp.cpp + src/GLATEdisp.cpp \ + src/GLATException.cpp # Define headers to be installed pkginclude_HEADERS = include/GLATObservation.hpp \ @@ -32,7 +33,8 @@ pkginclude_HEADERS = include/GLATObservation.hpp \ include/GLATRoi.hpp \ include/GLATResponse.hpp \ include/GLATResponseTable.hpp \ - include/GLATLib.hpp + include/GLATLib.hpp \ + include/GLATException.hpp # Build libtool library noinst_LTLIBRARIES = liblat.la diff --git a/inst/lat/include/GLATInstDir.hpp b/inst/lat/include/GLATInstDir.hpp index 351d8ef..fc9a5df 100644 --- a/inst/lat/include/GLATInstDir.hpp +++ b/inst/lat/include/GLATInstDir.hpp @@ -20,7 +20,7 @@ #define GLATINSTDIR_HPP /* __ Includes ___________________________________________________________ */ -#include +#include #include "GInstDir.hpp" #include "GSkyDir.hpp" @@ -35,46 +35,49 @@ ***************************************************************************/ class GLATInstDir : public GInstDir { - // I/O friends - friend std::ostream& operator<< (std::ostream& os, const GLATInstDir& dir); - public: // Constructors and destructors GLATInstDir(void); + explicit GLATInstDir(const GSkyDir& dir); GLATInstDir(const GLATInstDir& dir); virtual ~GLATInstDir(void); // Operators GLATInstDir& operator= (const GLATInstDir& dir); + // Implemented pure virtual base class methods + void clear(void); + GLATInstDir* clone(void) const; + std::string print(void) const; + // Methods - void clear(void); - void radec(const double& ra, const double& dec) { m_dir.radec(ra,dec); } - void radec_deg(const double& ra, const double& dec) { m_dir.radec_deg(ra,dec); } - void lb(const double& l, const double& b) { m_dir.lb(l,b); } - void lb_deg(const double& l, const double& b) { m_dir.lb_deg(l,b); } - double l(void) const { return m_dir.l(); } - double l_deg(void) const { return m_dir.l_deg(); } - double b(void) const { return m_dir.b(); } - double b_deg(void) const { return m_dir.b_deg(); } - double ra(void) const { return m_dir.ra(); } - double ra_deg(void) const { return m_dir.ra_deg(); } - double dec(void) const { return m_dir.dec(); } - double dec_deg(void) const { return m_dir.dec_deg(); } - double dist(GSkyDir& dir) const { return m_dir.dist(dir); } - double dist_deg(GSkyDir& dir) const { return m_dir.dist_deg(dir); } - double dist(GLATInstDir& dir) const; - double dist_deg(GLATInstDir& dir) const; + void skydir(const GSkyDir& dir) { m_dir=dir; } + void radec(const double& ra, const double& dec) { m_dir.radec(ra,dec); } + void radec_deg(const double& ra, const double& dec) { m_dir.radec_deg(ra,dec); } + void lb(const double& l, const double& b) { m_dir.lb(l,b); } + void lb_deg(const double& l, const double& b) { m_dir.lb_deg(l,b); } + GSkyDir skydir(void) const { return m_dir; } + double l(void) const { return m_dir.l(); } + double l_deg(void) const { return m_dir.l_deg(); } + double b(void) const { return m_dir.b(); } + double b_deg(void) const { return m_dir.b_deg(); } + double ra(void) const { return m_dir.ra(); } + double ra_deg(void) const { return m_dir.ra_deg(); } + double dec(void) const { return m_dir.dec(); } + double dec_deg(void) const { return m_dir.dec_deg(); } + double dist(GSkyDir& dir) const { return m_dir.dist(dir); } + double dist_deg(GSkyDir& dir) const { return m_dir.dist_deg(dir); } + double dist(GLATInstDir& dir) const; + double dist_deg(GLATInstDir& dir) const; protected: // Protected methods - void init_members(void); - void copy_members(const GLATInstDir& dir); - void free_members(void); - virtual GLATInstDir* clone(void) const; + void init_members(void); + void copy_members(const GLATInstDir& dir); + void free_members(void); // Data members - GSkyDir m_dir; //!< Incident direction of event + GSkyDir m_dir; //!< Observed incident direction of event }; #endif /* GLATINSTDIR_HPP */ diff --git a/inst/lat/pyext/latlib.i b/inst/lat/pyext/latlib.i index 20ae744..be0080b 100755 --- a/inst/lat/pyext/latlib.i +++ b/inst/lat/pyext/latlib.i @@ -17,7 +17,11 @@ /* __ LAT ________________________________________________________________ */ #include "GLATObservation.i" -//#include "GLATResponse.i" +#include "GLATEventCube.i" +#include "GLATEventBin.i" +#include "GLATResponse.i" +#include "GLATInstDir.i" +#include "GLATRoi.i" diff --git a/inst/lat/src/GLATInstDir.cpp b/inst/lat/src/GLATInstDir.cpp index b8bbb5c..812bbd7 100644 --- a/inst/lat/src/GLATInstDir.cpp +++ b/inst/lat/src/GLATInstDir.cpp @@ -20,8 +20,8 @@ #ifdef HAVE_CONFIG_H #include #endif -#include "GTools.hpp" #include "GLATInstDir.hpp" +#include "GTools.hpp" /* __ Method name definitions ____________________________________________ */ @@ -40,9 +40,9 @@ ==========================================================================*/ /***********************************************************************//** - * @brief Constructor + * @brief Void constructor ***************************************************************************/ -GLATInstDir::GLATInstDir(void) +GLATInstDir::GLATInstDir(void) : GInstDir() { // Initialise class members init_members(); @@ -53,12 +53,31 @@ GLATInstDir::GLATInstDir(void) /***********************************************************************//** + * @brief GSkyDir constructor + * + * @param[in] dir Sky direction. + * + * Construct LAT instrument direction from sky direction. + ***************************************************************************/ +GLATInstDir::GLATInstDir(const GSkyDir& dir) : GInstDir() +{ + // Initialise class members + init_members(); + + // Assign sky direction + m_dir = dir; + + // Return + return; +} + + +/***********************************************************************//** * @brief Copy constructor * - * @param[in] dir Instrument direction from which class should be - * instantiated. + * @param[in] dir Instrument direction. ***************************************************************************/ -GLATInstDir::GLATInstDir(const GLATInstDir& dir) +GLATInstDir::GLATInstDir(const GLATInstDir& dir) : GInstDir(dir) { // Initialise class members init_members(); @@ -93,13 +112,16 @@ GLATInstDir::~GLATInstDir(void) /***********************************************************************//** * @brief Assignment operator * - * @param[in] dir Instrument direction to be assigned. + * @param[in] dir Instrument direction. ***************************************************************************/ GLATInstDir& GLATInstDir::operator= (const GLATInstDir& dir) { // Execute only if object is not identical if (this != &dir) { + // Copy base class members + this->GInstDir::operator=(dir); + // Free members free_members(); @@ -123,14 +145,16 @@ GLATInstDir& GLATInstDir::operator= (const GLATInstDir& dir) ==========================================================================*/ /***********************************************************************//** - * @brief Clear sky direction + * @brief Clear instance ***************************************************************************/ void GLATInstDir::clear(void) { // Free members free_members(); + this->GInstDir::free_members(); // Initialise private members + this->GInstDir::init_members(); init_members(); // Return @@ -139,9 +163,18 @@ void GLATInstDir::clear(void) /***********************************************************************//** + * @brief Clone instance + ***************************************************************************/ +GLATInstDir* GLATInstDir::clone(void) const +{ + return new GLATInstDir(*this); +} + + +/***********************************************************************//** * @brief Compute angular distance between instrument directions in radians * - * @param[in] dir Instrument direction to which distance is to be computed. + * @param[in] dir Instrument direction. ***************************************************************************/ double GLATInstDir::dist(GLATInstDir& dir) const { @@ -162,7 +195,7 @@ double GLATInstDir::dist(GLATInstDir& dir) const /***********************************************************************//** * @brief Compute angular distance between instrument directions in degrees * - * @param[in] dir Instrument direction to which distance is to be computed. + * @param[in] dir Instrument direction. ***************************************************************************/ double GLATInstDir::dist_deg(GLATInstDir& dir) const { @@ -171,6 +204,22 @@ double GLATInstDir::dist_deg(GLATInstDir& dir) const } +/***********************************************************************//** + * @brief Print instrument direction information + ***************************************************************************/ +std::string GLATInstDir::print(void) const +{ + // Initialise result string + std::string result; + + // Append instrument direction + result.append("RA="+str(ra_deg())+", DEC="+str(dec_deg())); + + // Return result + return result; +} + + /*========================================================================== = = = Private methods = @@ -193,7 +242,7 @@ void GLATInstDir::init_members(void) /***********************************************************************//** * @brief Copy class members * - * @param[in] dir Instrument direction from which members should be copied + * @param[in] dir Instrument direction. ***************************************************************************/ void GLATInstDir::copy_members(const GLATInstDir& dir) { @@ -215,32 +264,8 @@ void GLATInstDir::free_members(void) } -/***********************************************************************//** - * @brief Clone instance - ***************************************************************************/ -GLATInstDir* GLATInstDir::clone(void) const -{ - return new GLATInstDir(*this); -} - - /*========================================================================== = = = Friends = = = ==========================================================================*/ - -/***********************************************************************//** - * @brief Output operator - * - * @param[in] os Output stream - * @param[in] column Instrument direction to put in output stream - ***************************************************************************/ -std::ostream& operator<< (std::ostream& os, const GLATInstDir& dir) -{ - // Output sky direction - os << dir.m_dir; - - // Return output stream - return os; -} diff --git a/inst/mwl/test/test_python.py b/inst/mwl/test/test_python.py index 1a62214..45e101b 100755 --- a/inst/mwl/test/test_python.py +++ b/inst/mwl/test/test_python.py @@ -87,6 +87,7 @@ def fit_spectrum(filename, xmlname): x = np.power(10., np.arange(-1., 8., 0.1)) y = [] d = GSkyDir() + d.radec_deg(83.6331, 22.0145) e = GEnergy() t = GTime() for energy in x: @@ -121,7 +122,8 @@ if __name__ == '__main__': extract_data("data/crab_mwl.fits") # Fit spectrum - fit_spectrum("data/crab_mwl.fits", "data/crab_mwl2.xml") +# fit_spectrum("data/crab_mwl.fits", "data/crab_mwl2.xml") + fit_spectrum("data/crab_mwl.fits", "data/crab_mwl3.xml") # Show spectrum plt.show() diff --git a/pyext/GOptimizerLM.i b/pyext/GOptimizerLM.i index f5eab3e..223cfeb 100644 --- a/pyext/GOptimizerLM.i +++ b/pyext/GOptimizerLM.i @@ -31,6 +31,7 @@ public: // Constructors and destructors GOptimizerLM(void); + GOptimizerLM(GLog& log); GOptimizerLM(const GOptimizerLM& opt); virtual ~GOptimizerLM(void); diff --git a/pyext/gammalib.i b/pyext/gammalib.i index 50a5055..f85d900 100755 --- a/pyext/gammalib.i +++ b/pyext/gammalib.i @@ -64,7 +64,7 @@ /* __ Application module _________________________________________________ */ #include "GApplication.i" -//#include "GLog.i" +#include "GLog.i" #include "GPars.i" #include "GPar.i" @@ -113,6 +113,7 @@ #include "GModelSpatialPtsrc.i" #include "GModelSpectral.i" #include "GModelSpectralPlaw.i" +#include "GModelSpectralExpPlaw.i" #include "GModelTemporal.i" #include "GModelTemporalConst.i" diff --git a/src/obs/GObservations_optimizer.cpp b/src/obs/GObservations_optimizer.cpp index b53004c..82a169e 100644 --- a/src/obs/GObservations_optimizer.cpp +++ b/src/obs/GObservations_optimizer.cpp @@ -34,8 +34,9 @@ /* __ Coding definitions _________________________________________________ */ /* __ Debug definitions __________________________________________________ */ -#define G_EVAL_TIMING 0 //!< Perform optimizer timing (0=no, 1=yes) -#define G_EVAL_DEBUG 0 //!< Perform optimizer debugging (0=no, 1=yes) +#define G_EVAL_TIMING 0 //!< Perform optimizer timing (0=no, 1=yes) +#define G_EVAL_DEBUG 0 //!< Perform optimizer debugging (0=no, 1=yes) +#define G_OPT_DEBUG 0 //!< Perform optimizer debugging (0=no, 1=yes) /* __ Prototypes _________________________________________________________ */ @@ -47,7 +48,7 @@ ==========================================================================*/ /***********************************************************************//** - * @brief Constructor + * @brief Void constructor ***************************************************************************/ GObservations::optimizer::optimizer(void) : GOptimizerFunction() { @@ -99,7 +100,7 @@ GObservations::optimizer::optimizer(const optimizer& fct) : GOptimizerFunction(f /***********************************************************************//** * @brief Destructor ***************************************************************************/ -GObservations::optimizer::~optimizer() +GObservations::optimizer::~optimizer(void) { // Free members free_members(); @@ -118,7 +119,7 @@ GObservations::optimizer::~optimizer() /***********************************************************************//** * @brief Assignment operator * - * @param[in] obs Instance to be assigned + * @param[in] obs Instance to be assigned. ***************************************************************************/ GObservations::optimizer& GObservations::optimizer::operator= (const optimizer& fct) { @@ -209,13 +210,17 @@ void GObservations::optimizer::eval(const GOptimizerPars& pars) m_npred += m_this->m_obs[i]->npred((GModels&)pars, m_wrk_grad); *m_gradient += *m_wrk_grad; #if G_EVAL_DEBUG - std::cout << "Npred=" << m_npred << " Grad=" - << *m_gradient << std::endl; + std::cout << "Unbinned Poisson"; + std::cout << " Npred=" << m_npred; + std::cout << " Grad="<< *m_gradient << std::endl; #endif // Update the log-likelihood poisson_unbinned(*(m_this->m_obs[i]), pars); + // Add the Npred value to the log-likelihood + m_value += m_npred; + } // endif: Poisson statistics // ... otherwise throw an exception @@ -229,12 +234,20 @@ void GObservations::optimizer::eval(const GOptimizerPars& pars) else { // Poisson statistics - if (toupper(statistics) == "POISSON") + if (toupper(statistics) == "POISSON") { + #if G_EVAL_DEBUG + std::cout << "Binned Poisson" << std::endl; + #endif poisson_binned(*(m_this->m_obs[i]), pars); + } // ... or Gaussian statistics - else if (toupper(statistics) == "GAUSSIAN") + else if (toupper(statistics) == "GAUSSIAN") { + #if G_EVAL_DEBUG + std::cout << "Binned Gaussian" << std::endl; + #endif gaussian_binned(*(m_this->m_obs[i]), pars); + } // ... or unsupported else @@ -245,9 +258,6 @@ void GObservations::optimizer::eval(const GOptimizerPars& pars) } // endfor: looped over observations - // Add the Npred value to the log-likelihood - m_value += m_npred; - // Release stack m_covar->stack_destroy(); @@ -311,9 +321,7 @@ void GObservations::optimizer::poisson_unbinned(const GObservation& obs, GEvent* event = obs.events()->pointer(i); // Get model and derivative - double model = obs.model((GModels&)pars, - event->dir(), event->energy(), event->time(), - m_wrk_grad); + double model = obs.model((GModels&)pars, *event, m_wrk_grad); // Skip bin if model is too small (avoids -Inf or NaN gradients) if (model <= m_minmod) @@ -410,6 +418,17 @@ void GObservations::optimizer::poisson_binned(const GObservation& obs, clock_t t_start = clock(); #endif + // Initialise statistics + #if G_OPT_DEBUG + int n_bins = 0; + int n_used = 0; + int n_small_model = 0; + int n_zero_data = 0; + double sum_data = 0.0; + double sum_model = 0.0; + double init_value = m_value; + #endif + // Get number of parameters int npars = pars.npars(); @@ -420,6 +439,11 @@ void GObservations::optimizer::poisson_binned(const GObservation& obs, // Iterate over all bins for (int i = 0; i < obs.events()->size(); ++i) { + // Update number of bins + #if G_OPT_DEBUG + n_bins++; + #endif + // Get pointer to bin GEventBin* bin = (GEventBin*)obs.events()->pointer(i); @@ -427,16 +451,25 @@ void GObservations::optimizer::poisson_binned(const GObservation& obs, double data = bin->counts(); // Get model and derivative - double model = obs.model((GModels&)pars, - bin->dir(), bin->energy(), bin->time(), - m_wrk_grad); + double model = obs.model((GModels&)pars, *bin, m_wrk_grad); // Multiply model by bin size model *= bin->size(); // Skip bin if model is too small (avoids -Inf or NaN gradients) - if (model <= m_minmod) + if (model <= m_minmod) { + #if G_OPT_DEBUG + n_small_model++; + #endif continue; + } + + // Update statistics + #if G_OPT_DEBUG + n_used++; + sum_data += data; + sum_model += model; + #endif // Update Npred m_npred += model; @@ -498,6 +531,11 @@ void GObservations::optimizer::poisson_binned(const GObservation& obs, // ... handle now data=0 else { + // Update statistics + #if G_OPT_DEBUG + n_zero_data++; + #endif + // Update Poissonian statistics (excluding factorial term for // faster computation) m_value += model; @@ -519,6 +557,18 @@ void GObservations::optimizer::poisson_binned(const GObservation& obs, if (values != NULL) delete [] values; if (inx != NULL) delete [] inx; + // Dump statistics + #if G_OPT_DEBUG + std::cout << "Number of bins: " << n_bins << std::endl; + std::cout << "Number of bins used for computation: " << n_used << std::endl; + std::cout << "Number of bins excluded due to small model: " << n_small_model << std::endl; + std::cout << "Number of bins with zero data: " << n_zero_data << std::endl; + std::cout << "Sum of data: " << sum_data << std::endl; + std::cout << "Sum of model: " << sum_model << std::endl; + std::cout << "Initial statistics: " << init_value << std::endl; + std::cout << "Statistics: " << m_value-init_value << std::endl; + #endif + // Optionally dump gradient and covariance matrix #if G_EVAL_DEBUG std::cout << *m_gradient << std::endl; @@ -589,9 +639,7 @@ void GObservations::optimizer::gaussian_binned(const GObservation& obs, continue; // Get model and derivative - double model = obs.model((GModels&)pars, - bin->dir(), bin->energy(), bin->time(), - m_wrk_grad); + double model = obs.model((GModels&)pars, *bin, m_wrk_grad); // Multiply model by bin size model *= bin->size();