Simon Bonnefoy / gammalib

Commit e7028dbbd1e5bb0ffbd2dd694d5c89a4d09713e1

Authored by Jürgen Knödlseder
1 parent f0517f5c

Integrate code

Inline Side-by-side
Showing 15 changed files with 563 additions and 508 deletions
ChangeLog
  View file @e7028db
1   -2017-12-01
  1 +2018-01-09
2 2  
3 3 * Version 1.5.0 released
4 4 =========================
5 5  
  6 + Add GPha::counts_spectrum() and GPha::backscal_spectrum() methods (#2290)
  7 + Add GNdarray sign(), operator*= and operator/= (#2290)
  8 + Add GSkyMap::counts() and GSkyMap::flux() methods and clip() friend operator (#2290)
6 9 Add GSkyMap::array() method to Python interface (#2285)
7 10 CTA observations handle "statistic" attribute in XML file (#2238)
8 11 Use CSTAT and CHI2 instead of POISSON and GAUSSIAN for statistic (#2237)
... ...
1 1 New Features and Important Changes in GammaLib 1.5.0
2 2  
3   -1 December 2017
  3 +9 January 2018
4 4  
5 5  
6 6 1. Introduction
... ... @@ -45,6 +45,9 @@ The following methods have been added:
45 45 - GEbounds::operator!=(GEbounds&, GEbounds&)
46 46 - GSkyMap::smooth(std::string&, double&)
47 47 - GSkyMap::array(int&) - only for Python interface
  48 +- GSkyMap::counts()
  49 +- GSkyMap::flux()
  50 +- clip(GSkyMap&, double&)
48 51 - GNdarray::shape(std::vector<int>&)
49 52 - GPha::operator+=(GPha&)
50 53 - GPha::operator-=(GPha&)
... ... @@ -54,6 +57,8 @@ The following methods have been added:
54 57 - GPha::areascal(int&)
55 58 - GPha::backscal(int&, double&)
56 59 - GPha::backscal(int&)
  60 +- GPha::counts_spectrum()
  61 +- GPha::backscal_spectrum()
57 62 - operator+(GPha&, GPha&)
58 63 - operator-(GPha&, GPha&)
59 64 - operator*(GPha&, double&)
... ... @@ -135,7 +140,8 @@ conversion of any spectral model into a spectral node model (#2178).
135 140  
136 141 9. Numerics module
137 142 ------------------
138   -None
  143 +Added GNdarray sign() friend operator and operator*= and operator/=
  144 +(#2290).
139 145  
140 146  
141 147 10. Observation module
... ... @@ -166,6 +172,9 @@ None
166 172  
167 173 12. Sky module
168 174 --------------
  175 +Added GSkyMap::counts() and GSkyMap::flux() methods and clip() frield
  176 +operators that allows clipping of sky map values (#2290)
  177 +
169 178 Added the GSkyMap::array() method to the Python interface to returns the sky
170 179 map pixels as a list of lists (#2285)
171 180  
... ... @@ -196,6 +205,9 @@ None
196 205  
197 206 16. Xspec module
198 207 ----------------
  208 +Added GPha::counts_spectrum() and GPha::backscal_spectrum() methods to access
  209 +the GPha content (#2290).
  210 +
199 211 Added arithmetic operators for GPha class to add, subtract, scale and divide
200 212 spectra, and methods to set and retrieve the area and background scaling
201 213 factors.
... ...
doc/source/admin/release_history/1.5.rst
  View file @e7028db
... ... @@ -30,6 +30,10 @@ Bug fixes
30 30 Improvements
31 31 ------------
32 32  
  33 +* [`2290 <https://cta-redmine.irap.omp.eu/issues/2290>`_] -
  34 + Add ``GSkyMap::counts()`` and ``GSkyMap::flux()`` methods and ``clip()``
  35 + friend operator, ``GNdarray`` ``sign()``, ``operator*=`` and ``operator/=``,
  36 + and ``GPha::counts_spectrum()`` and ``GPha::backscal_spectrum()`` methods
33 37 * [`2285 <https://cta-redmine.irap.omp.eu/issues/2285>`_] -
34 38 Add ``GSkyMap::array()`` method to Python interface
35 39 * [`2238 <https://cta-redmine.irap.omp.eu/issues/2238>`_] -
... ...
include/GNdarray.hpp
  View file @e7028db
1 1 /***************************************************************************
2 2 * GNdarray.hpp - N-dimensional array class *
3 3 * ----------------------------------------------------------------------- *
4   - * copyright (C) 2016-2017 by Juergen Knoedlseder *
  4 + * copyright (C) 2016-2018 by Juergen Knoedlseder *
5 5 * ----------------------------------------------------------------------- *
6 6 * *
7 7 * This program is free software: you can redistribute it and/or modify *
... ...
include/GPha.hpp
  View file @e7028db
1 1 /***************************************************************************
2 2 * GPha.hpp - XSPEC Pulse Height Analyzer class *
3 3 * ----------------------------------------------------------------------- *
4   - * copyright (C) 2013-2017 by Juergen Knoedlseder *
  4 + * copyright (C) 2013-2018 by Juergen Knoedlseder *
5 5 * ----------------------------------------------------------------------- *
6 6 * *
7 7 * This program is free software: you can redistribute it and/or modify *
... ... @@ -122,7 +122,7 @@ protected:
122 122 void copy_members(const GPha& pha);
123 123 void free_members(void);
124 124 void alloc(const int& size);
125   -
  125 +
126 126 // Protected members
127 127 mutable GFilename m_filename; //!< Filename of origin
128 128 std::vector<double> m_counts; //!< Counts data
... ...
include/GSkyMap.hpp
  View file @e7028db
1 1 /***************************************************************************
2 2 * GSkyMap.hpp - Sky map class *
3 3 * ----------------------------------------------------------------------- *
4   - * copyright (C) 2010-2017 by Juergen Knoedlseder *
  4 + * copyright (C) 2010-2018 by Juergen Knoedlseder *
5 5 * ----------------------------------------------------------------------- *
6 6 * *
7 7 * This program is free software: you can redistribute it and/or modify *
... ... @@ -92,7 +92,7 @@ class GSkyMap : public GBase {
92 92 friend GSkyMap log10(const GSkyMap& map);
93 93 friend GSkyMap abs(const GSkyMap& map);
94 94 friend GSkyMap sign(const GSkyMap& map);
95   - friend GSkyMap clip(const GSkyMap& map,const double& thresh);
  95 + friend GSkyMap clip(const GSkyMap& map, const double& thresh);
96 96  
97 97 public:
98 98 // Constructors and destructors
... ... @@ -158,8 +158,10 @@ public:
158 158 int pix2inx(const GSkyPixel& pixel) const;
159 159 int dir2inx(const GSkyDir& dir) const;
160 160 GSkyPixel dir2pix(const GSkyDir& dir) const;
  161 + GNdarray counts(void) const;
161 162 double flux(const int& index, const int& map = 0) const;
162 163 double flux(const GSkyPixel& pixel, const int& map = 0) const;
  164 + GNdarray flux(void) const;
163 165 double solidangle(const int& index) const;
164 166 double solidangle(const GSkyPixel& pixel) const;
165 167 bool contains(const GSkyDir& dir) const;
... ... @@ -171,8 +173,6 @@ public:
171 173 const double* pixels(void) const;
172 174 GSkyMap extract(const int& map, const int& nmaps = 1) const;
173 175 void stack_maps(void);
174   - GNdarray total_counts(void) const;
175   - GNdarray total_flux(void) const;
176 176 void load(const GFilename& filename);
177 177 void save(const GFilename& filename,
178 178 const bool& clobber = false) const;
... ...
inst/cta/include/GCTAOnOffObservation.hpp
  View file @e7028db
1 1 /***************************************************************************
2 2 * GCTAOnOffObservation.hpp - CTA On/Off observation class *
3 3 * ----------------------------------------------------------------------- *
4   - * copyright (C) 2013-2017 by Chia-Chun Lu & Christoph Deil *
  4 + * copyright (C) 2013-2018 by Chia-Chun Lu & Christoph Deil *
5 5 * ----------------------------------------------------------------------- *
6 6 * *
7 7 * This program is free software: you can redistribute it and/or modify *
... ... @@ -140,14 +140,6 @@ protected:
140 140 double N_bgd(const GModels& models, const int& ibin, GVector* grad) const;
141 141  
142 142 // Likelihood methods
143   - double wstat_value(double non,
144   - double noff,
145   - double alpha,
146   - double ngam,
147   - double& nonpred,
148   - double& nbgd,
149   - double& dlogLdsky,
150   - double& d2logLdsky2) const;
151 143 virtual double likelihood_cstat(const GModels& models,
152 144 GVector* gradient,
153 145 GMatrixSparse* curvature,
... ... @@ -156,6 +148,14 @@ protected:
156 148 GVector* gradient,
157 149 GMatrixSparse* curvature,
158 150 double* npred) const;
  151 + virtual double wstat_value(const double& non,
  152 + const double& noff,
  153 + const double& alpha,
  154 + const double& ngam,
  155 + double& nonpred,
  156 + double& nbgd,
  157 + double& dlogLdsky,
  158 + double& d2logLdsky2) const;
159 159  
160 160 // Energy dispersion integration kernel
161 161 class edisp_kern : public GFunction {
... ...
inst/cta/src/GCTAOnOffObservation.cpp
  View file @e7028db
1 1 /***************************************************************************
2 2 * GCTAOnOffObservation.cpp - CTA On/Off observation class *
3 3 * ----------------------------------------------------------------------- *
4   - * copyright (C) 2013-2017 by Chia-Chun Lu & Christoph Deil *
  4 + * copyright (C) 2013-2018 by Chia-Chun Lu & Christoph Deil *
5 5 * ----------------------------------------------------------------------- *
6 6 * *
7 7 * This program is free software: you can redistribute it and/or modify *
... ... @@ -70,12 +70,6 @@ const GObservationRegistry g_onoff_obs_cta_registry(&amp;g_onoff_obs_cta_seed);
70 70 #define G_LIKELIHOOD "GCTAOnOffObservation::likelihood(GModels&, "\
71 71 "GOptimizerPars&, GMatrixSparse&, "\
72 72 "GVector&, double&, double&)"
73   -#define G_LIKELIHOOD_CSTAT "GCTAOnOffObservation::likelihood_cstat("\
74   - "GModels&, GOptimizerPars&, GMatrixSparse&, "\
75   - "GVector&, double&, double&)"
76   -#define G_LIKELIHOOD_WSTAT "GCTAOnOffObservation::likelihood_wstat("\
77   - "GModels&, GOptimizerPars&, GMatrixSparse&, "\
78   - "GVector&, double&, double&)"
79 73 #define G_SET "GCTAOnOffObservation::set(GCTAObservation&, GSkyDir&, "\
80 74 "GSkyRegionMap&, GSkyRegionMap&)"
81 75 #define G_COMPUTE_ARF "GCTAOnOffObservation::compute_arf(GCTAObservation&, "\
... ... @@ -86,6 +80,13 @@ const GObservationRegistry g_onoff_obs_cta_registry(&amp;g_onoff_obs_cta_seed);
86 80 "GCTAObservation&, GSkyRegionMap&, GSkyRegionMap&)"
87 81 #define G_COMPUTE_RMF "GCTAOnOffObservation::compute_rmf(GCTAObservation&, "\
88 82 "GSkyRegionMap&)"
  83 +#define G_MODEL_BACKGROUND "GCTAOnOffObservation::model_background(GModels&)"
  84 +#define G_LIKELIHOOD_CSTAT "GCTAOnOffObservation::likelihood_cstat("\
  85 + "GModels&, GOptimizerPars&, GMatrixSparse&, "\
  86 + "GVector&, double&, double&)"
  87 +#define G_LIKELIHOOD_WSTAT "GCTAOnOffObservation::likelihood_wstat("\
  88 + "GModels&, GOptimizerPars&, GMatrixSparse&, "\
  89 + "GVector&, double&, double&)"
89 90  
90 91 /* __ Constants __________________________________________________________ */
91 92 const double minmod = 1.0e-100; //!< Minimum model value
... ... @@ -1736,36 +1737,36 @@ double GCTAOnOffObservation::N_bgd(const GModels&amp; models,
1736 1737 }
1737 1738  
1738 1739 /***********************************************************************
1739   - * @brief Compute predicted gamma-ray counts Pha for given model
  1740 + * @brief Compute predicted gamma-ray counts for given model
1740 1741 *
1741 1742 * @param[in] models Model container.
1742   - * @return model gamma Pha.
1743   - *
1744   - * Returns spectrum of predicted number of source events in the On regions
  1743 + * @return model Pha with predicted gamma-ray counts.
1745 1744 *
  1745 + * Returns spectrum of predicted number of source events in the On
  1746 + * regions.
1746 1747 ***********************************************************************/
1747 1748 GPha GCTAOnOffObservation::model_gamma(const GModels& models) const
1748 1749 {
1749   - // Get number of model parameters in model container
1750   - int npars = models.npars();
  1750 + // Get number of model parameters in model container
  1751 + int npars = models.npars();
1751 1752  
1752   - // Create dummy gradient vectors to provide required argument to N_gamma
1753   - GVector dummy_grad(npars);
  1753 + // Create dummy gradient vectors to provide required argument to N_gamma
  1754 + GVector dummy_grad(npars);
1754 1755  
1755   - // Initialise output Pha
1756   - GEbounds ereco = m_on_spec.ebounds();
1757   - GPha gammas = GPha(ereco);
  1756 + // Initialise output spectrum
  1757 + GEbounds ereco = m_on_spec.ebounds();
  1758 + GPha gammas = GPha(ereco);
1758 1759  
1759   - // Loop over all energy bins
1760   - for (int i = 0; i < m_on_spec.size(); ++i) {
1761   - double ngam = N_gamma(models, i, &dummy_grad);
1762   - gammas.fill(ereco.emean(i),ngam);
1763   - }
  1760 + // Loop over all energy bins
  1761 + for (int i = 0; i < m_on_spec.size(); ++i) {
  1762 + gammas.fill(ereco.emean(i), N_gamma(models, i, &dummy_grad));
  1763 + }
1764 1764  
1765   - // Return model coun pha
1766   - return gammas;
  1765 + // Return predicted gamma-ray counts spectrum
  1766 + return gammas;
1767 1767 }
1768 1768  
  1769 +
1769 1770 /***********************************************************************
1770 1771 * @brief Compute predicted background counts Pha for given model
1771 1772 *
... ... @@ -1781,70 +1782,73 @@ GPha GCTAOnOffObservation::model_gamma(const GModels&amp; models) const
1781 1782 ***********************************************************************/
1782 1783 GPha GCTAOnOffObservation::model_background(const GModels& models) const
1783 1784 {
1784   - // Get number of model parameters in model container
1785   - int npars = models.npars();
  1785 + // Get number of model parameters in model container
  1786 + int npars = models.npars();
1786 1787  
1787   - // Create dummy gradient vectors to provide required argument to N_gamma
1788   - GVector dummy_grad(npars);
  1788 + // Create dummy gradient vectors to provide required argument to N_gamma
  1789 + GVector dummy_grad(npars);
1789 1790  
1790   - // Initialise output Pha
1791   - GEbounds ereco = m_on_spec.ebounds();
1792   - GPha bgds = GPha(ereco);
  1791 + // Initialise output spectrum
  1792 + GEbounds ereco = m_on_spec.ebounds();
  1793 + GPha bgds = GPha(ereco);
1793 1794  
1794   - // Extract statistic for this observation
1795   - std::string statistic = gammalib::toupper(this->statistic());
  1795 + // Extract statistic for this observation
  1796 + std::string statistic = gammalib::toupper(this->statistic());
1796 1797  
1797   - // Loop over all energy bins
1798   - for (int i = 0; i < m_on_spec.size(); ++i) {
  1798 + // Loop over all energy bins
  1799 + for (int i = 0; i < m_on_spec.size(); ++i) {
1799 1800  
1800   - // Initialise variable to store number of background counts
1801   - double nbgd = 0.0;
  1801 + // Initialise variable to store number of background counts
  1802 + double nbgd = 0.0;
1802 1803  
1803   - // Choose background evaluation method based on fit statistics
  1804 + // Choose background evaluation method based on fit statistics
1804 1805  
1805   - // CSTAT
1806   - if ((statistic == "POISSON") || (statistic == "CSTAT")) {
1807   - nbgd = N_bgd(models, i, &dummy_grad);
1808   - }
1809   -
1810   - // ... or Poisson statistic with measured background
1811   - else if (statistic == "WSTAT") {
1812   - // Get number of On and Off counts
1813   - double non = m_on_spec[i];
1814   - double noff = m_off_spec[i];
1815   -
1816   - // Get background scaling
1817   - double alpha = m_on_spec.backscal(i);
1818   -
1819   - // Get number of gamma and background events (and corresponding
1820   - // spectral model gradients)
1821   - double ngam = N_gamma(models, i, &dummy_grad);
1822   -
1823   - // Initialise variables for likelihood computation
1824   - double nonpred = 0.0;
1825   - double dlogLdsky = 0.0;
1826   - double d2logLdsky2 = 0.0;
1827   -
1828   - // Perform likelihood profiling and derive likelihood value
1829   - // The number of background counts is updated to the profile value
1830   - double logL = wstat_value(non,noff,alpha,ngam,
1831   - nonpred,nbgd,dlogLdsky,d2logLdsky2);
1832   - }
1833   -
1834   - // ... or unsupported
1835   - else {
1836   - std::string msg = "Invalid statistic \""+statistic+"\" encountered. "
1837   - "Either specify \"POISSON\", \"CSTAT\" or "
1838   - "\"WSTAT\".";
1839   - throw GException::invalid_value(G_LIKELIHOOD, msg);
1840   - }
1841   -
1842   - // Fill background Pha
1843   - bgds.fill(ereco.emean(i),nbgd);
1844   - }
1845   -
1846   - // Return model coun pha
1847   - return bgds;
  1806 + // CSTAT
  1807 + if ((statistic == "POISSON") || (statistic == "CSTAT")) {
  1808 + nbgd = N_bgd(models, i, &dummy_grad);
  1809 + }
  1810 +
  1811 + // ... or Poisson statistic with measured background
  1812 + else if (statistic == "WSTAT") {
  1813 +
  1814 + // Get number of On and Off counts
  1815 + double non = m_on_spec[i];
  1816 + double noff = m_off_spec[i];
  1817 +
  1818 + // Get background scaling
  1819 + double alpha = m_on_spec.backscal(i);
  1820 +
  1821 + // Get number of gamma and background events (and corresponding
  1822 + // spectral model gradients)
  1823 + double ngam = N_gamma(models, i, &dummy_grad);
  1824 +
  1825 + // Initialise variables for likelihood computation
  1826 + double nonpred = 0.0;
  1827 + double dlogLdsky = 0.0;
  1828 + double d2logLdsky2 = 0.0;
  1829 +
  1830 + // Perform likelihood profiling and derive likelihood value
  1831 + // The number of background counts is updated to the profile value
  1832 + double logL = wstat_value(non, noff, alpha, ngam, nonpred, nbgd,
  1833 + dlogLdsky,d2logLdsky2);
  1834 +
  1835 + } // endelse: WSTAT
  1836 +
  1837 + // ... or unsupported
  1838 + else {
  1839 + std::string msg = "Invalid statistic \""+statistic+"\" encountered. "
  1840 + "Either specify \"POISSON\", \"CSTAT\" or "
  1841 + "\"WSTAT\".";
  1842 + throw GException::invalid_value(G_MODEL_BACKGROUND, msg);
  1843 + }
  1844 +
  1845 + // Fill background spectrum
  1846 + bgds.fill(ereco.emean(i), nbgd);
  1847 +
  1848 + } // endfor: looped over energy bins
  1849 +
  1850 + // Return model count spectrum
  1851 + return bgds;
1848 1852 }
1849 1853  
1850 1854  
... ... @@ -1951,9 +1955,9 @@ GPha GCTAOnOffObservation::model_background(const GModels&amp; models) const
1951 1955 * \f]
1952 1956 ***********************************************************************/
1953 1957 double GCTAOnOffObservation::likelihood_cstat(const GModels& models,
1954   - GVector* gradient,
1955   - GMatrixSparse* curvature,
1956   - double* npred) const
  1958 + GVector* gradient,
  1959 + GMatrixSparse* curvature,
  1960 + double* npred) const
1957 1961 {
1958 1962 // Timing measurement
1959 1963 #if defined(G_LIKELIHOOD_DEBUG)
... ... @@ -2174,6 +2178,188 @@ double GCTAOnOffObservation::likelihood_cstat(const GModels&amp; models,
2174 2178 return value;
2175 2179 }
2176 2180  
  2181 +
  2182 +
  2183 +/***********************************************************************//**
  2184 + * @brief Evaluate log-likelihood function for On/Off analysis in the
  2185 + * case of Poisson signal with measured Poisson background
  2186 + *
  2187 + * @param[in] models Models.
  2188 + * @param[in,out] gradient Pointer to gradients.
  2189 + * @param[in,out] curvature Pointer to curvature matrix.
  2190 + * @param[in,out] npred Pointer to Npred value.
  2191 + * @return Log-likelihood value.
  2192 + *
  2193 + * @exception GException::invalid_value
  2194 + * There are no model parameters.
  2195 + *
  2196 + * Computes the log-likelihood value for the On/Off observation. The method
  2197 + * loops over the energy bins to update the function value, its derivatives
  2198 + * and the curvature matrix.
  2199 + *
  2200 + * The number of On counts \f$n_{\rm on}\f$ and Off counts \f$n_{\rm off}\f$
  2201 + * are taken from the On and Off spectra, the expected number of gamma-ray
  2202 + * events \f$\mu_s\f$ is computed from the spectral models of the relevant
  2203 + * components in the model container (spatial and temporal components
  2204 + * are ignored so far). See the N_gamma() method for details about the
  2205 + * model computations.
  2206 + *
  2207 + * The estimated number of background counts \f$\mu_b\f$ is derived based
  2208 + * on the measurement in the Off region by analytical minimization of
  2209 + * the Poisson likelihood, i.e., it is treated as a nuisance parameter.
  2210 + * See the wstat_value function for the derivation.
  2211 + ***********************************************************************/
  2212 +double GCTAOnOffObservation::likelihood_wstat(const GModels& models,
  2213 + GVector* gradient,
  2214 + GMatrixSparse* curvature,
  2215 + double* npred) const
  2216 +{
  2217 + // Timing measurement
  2218 + #if defined(G_LIKELIHOOD_DEBUG)
  2219 + #ifdef _OPENMP
  2220 + double t_start = omp_get_wtime();
  2221 + #else
  2222 + clock_t t_start = clock();
  2223 + #endif
  2224 + #endif
  2225 +
  2226 + // Initialise statistics
  2227 + #if defined(G_LIKELIHOOD_DEBUG)
  2228 + int n_bins = m_on_spec.size();
  2229 + int n_used = 0;
  2230 + int n_small_model = 0;
  2231 + int n_zero_data = 0;
  2232 + double sum_data = 0.0;
  2233 + double sum_model = 0.0;
  2234 + double init_npred = *npred;
  2235 + #endif
  2236 +
  2237 + // Initialise log-likelihood value
  2238 + double value = 0.0;
  2239 +
  2240 + // Get number of model parameters in model container
  2241 + int npars = models.npars();
  2242 +
  2243 + // Create model gradient vectors for sky parameters
  2244 + GVector sky_grad(npars);
  2245 +
  2246 + // Allocate working array
  2247 + GVector colvar(npars);
  2248 +
  2249 + // Check that there is at least one parameter
  2250 + if (npars > 0) {
  2251 +
  2252 + // Loop over all energy bins
  2253 + for (int i = 0; i < m_on_spec.size(); ++i) {
  2254 +
  2255 + // Reinitialize working arrays
  2256 + for (int j = 0; j < npars; ++j) {
  2257 + sky_grad[j] = 0.0;
  2258 + }
  2259 +
  2260 + // Get number of On and Off counts
  2261 + double non = m_on_spec[i];
  2262 + double noff = m_off_spec[i];
  2263 +
  2264 + // Get background scaling
  2265 + double alpha = m_on_spec.backscal(i);
  2266 +
  2267 + // Get number of gamma and background events (and corresponding
  2268 + // spectral model gradients)
  2269 + double ngam = N_gamma(models, i, &sky_grad);
  2270 +
  2271 + // Initialise variables for likelihood computation
  2272 + double nonpred = 0.0;
  2273 + double nbgd = 0.0;
  2274 + double dlogLdsky = 0.0;
  2275 + double d2logLdsky2 = 0.0;
  2276 +
  2277 + // Perform likelihood profiling and derive likelihood value
  2278 + double logL = wstat_value(non, noff, alpha, ngam,
  2279 + nonpred, nbgd, dlogLdsky, d2logLdsky2);
  2280 +
  2281 + // Update global log-likelihood and Npred
  2282 + value += logL;
  2283 + *npred += nonpred;
  2284 +
  2285 + // Update statistics
  2286 + #if defined(G_LIKELIHOOD_DEBUG)
  2287 + n_used++;
  2288 + sum_data += non;
  2289 + sum_model += nonpred;
  2290 + #endif
  2291 +
  2292 + // Loop over all parameters
  2293 + for (int j = 0; j < npars; ++j) {
  2294 +
  2295 + // If spectral model for sky component is non-zero and
  2296 + // non-infinite then handle sky component gradients and
  2297 + // second derivatives including at least a sky component ...
  2298 + if (sky_grad[j] != 0.0 && !gammalib::is_infinite(sky_grad[j])) {
  2299 +
  2300 + // Gradient
  2301 + (*gradient)[j] += dlogLdsky * sky_grad[j];
  2302 +
  2303 + // Hessian (from first-order derivatives only)
  2304 + for (int k = 0; k < npars; ++k) {
  2305 +
  2306 + // If spectral model for sky component is non-zero and
  2307 + // non-infinite then we have the curvature element
  2308 + // of a sky component
  2309 + if (sky_grad[k] != 0.0 &&
  2310 + !gammalib::is_infinite(sky_grad[k])) {
  2311 + colvar[k] = sky_grad[j] * sky_grad[k] * d2logLdsky2;
  2312 + }
  2313 +
  2314 + // ...... else if spectral model for background component
  2315 + // or neither sky nor background
  2316 + else {
  2317 + colvar[k] = 0.0;
  2318 + }
  2319 +
  2320 + } // endfor: Hessian computation
  2321 +
  2322 + // Update matrix
  2323 + curvature->add_to_column(j, colvar);
  2324 +
  2325 + } // endif: spectral model is non-zero and non-infinite
  2326 +
  2327 + } // endfor: looped over all parameters for derivatives computation
  2328 +
  2329 + } // endfor: looped over energy bins
  2330 +
  2331 + } // endif: number of parameters was positive
  2332 +
  2333 + // ... else there are no parameters, so throw an exception
  2334 + else {
  2335 + std::string msg ="No model parameter for the computation of the "
  2336 + "likelihood in observation \""+this->name()+
  2337 + "\" (ID "+this->id()+").\n";
  2338 + throw GException::invalid_value(G_LIKELIHOOD_WSTAT, msg);
  2339 + }
  2340 +
  2341 + // Optionally dump gradient and curvature matrix
  2342 + #if defined(G_LIKELIHOOD_DEBUG)
  2343 + std::cout << *gradient << std::endl;
  2344 + std::cout << *curvature << std::endl;
  2345 + #endif
  2346 +
  2347 + // Timing measurement
  2348 + #if defined(G_LIKELIHOOD_DEBUG)
  2349 + #ifdef _OPENMP
  2350 + double t_elapse = omp_get_wtime()-t_start;
  2351 + #else
  2352 + double t_elapse = (double)(clock() - t_start) / (double)CLOCKS_PER_SEC;
  2353 + #endif
  2354 + std::cout << "GCTAOnOffObservation::optimizer::likelihood_wstat: CPU usage = "
  2355 + << t_elapse << " sec" << std::endl;
  2356 + #endif
  2357 +
  2358 + // Return
  2359 + return value;
  2360 +}
  2361 +
  2362 +
2177 2363 /***********************************************************************//**
2178 2364 * @brief Evaluate log-likelihood value in energy bin for On/Off analysis
2179 2365 * by profiling over the number of background counts
... ... @@ -2379,279 +2565,96 @@ double GCTAOnOffObservation::likelihood_cstat(const GModels&amp; models,
2379 2565 * overestimated if for some bins \f$n_{\rm on}=0\f$ and/or
2380 2566 * \f$n_{\rm off}=0\f$ (i.e. if the special cases are encountered).
2381 2567 ***********************************************************************/
2382   -
2383   -double GCTAOnOffObservation::wstat_value(double non,
2384   - double noff,
2385   - double alpha,
2386   - double ngam,
2387   - double& nonpred,
2388   - double& nbgd,
2389   - double& dlogLdsky,
2390   - double& d2logLdsky2) const
  2568 +double GCTAOnOffObservation::wstat_value(const double& non,
  2569 + const double& noff,
  2570 + const double& alpha,
  2571 + const double& ngam,
  2572 + double& nonpred,
  2573 + double& nbgd,
  2574 + double& dlogLdsky,
  2575 + double& d2logLdsky2) const
2391 2576 {
2392   - // Initialise log-likelihood value
2393   - double logL;
2394   -
2395   - // Precompute some values
2396   - double alphap1 = alpha + 1.0;
2397   - double alpharat = alpha / alphap1;
2398   -
2399   - // Calculate number of background events, profile likelihood value
2400   - // and likelihood derivatives
2401   -
2402   - // Special case noff = 0
2403   - if (noff == 0.0) {
2404   -
2405   - // Case A: non = 0. In this case nbgd < 0 hence we set nbgd = 0
2406   - if (non == 0.0) {
2407   - nbgd = 0;
2408   - nonpred = ngam;
2409   - logL = ngam;
2410   - dlogLdsky = 1.0;
2411   - d2logLdsky2 = 0.0;
2412   - }
2413   -
2414   - // Case B: nbgd is positive
2415   - else if (ngam < non * alpharat) {
2416   - nbgd = non / alphap1 - ngam / alpha;
2417   - nonpred = ngam + alpha * nbgd;
2418   - logL = -ngam / alpha - non * std::log(alpharat);
2419   - dlogLdsky = -1.0/alpha;
2420   - d2logLdsky2 = 0.0;
2421   - }
2422   -
2423   - // Case C: nbgd is zero or negative, hence set nbgd = 0
2424   - else {
2425   - nbgd = 0;
2426   - nonpred = ngam;
2427   - logL = ngam + non * (std::log(non) - std::log(ngam) - 1.0);
2428   - dlogLdsky = 1.0 - non / ngam;
2429   - d2logLdsky2 = non / (ngam * ngam);
2430   - }
2431   -
2432   - } // endif: noff = 0
2433   -
2434   - // Special case non = 0
2435   - else if (non == 0.0) {
2436   - nbgd = 0;
2437   - nonpred = ngam + alpharat * noff;
2438   - logL = ngam + noff * std::log(alphap1);
2439   - dlogLdsky = 1.0;
2440   - d2logLdsky2 = 0.0;
2441   - } // endif: non = 0
2442   -
2443   - // General case
2444   - else {
2445   -
2446   - // Compute log-likelihood value
2447   - double alphat2 = 2.0 * alpha;
2448   - double C = alpha * (non + noff) - alphap1 * ngam;
2449   - double D = std::sqrt(C*C + 4.0 * alpha * alphap1 * noff * ngam);
2450   - nbgd = (C + D) / (alphat2 * alphap1);
2451   - nonpred = ngam + alpha * nbgd;
2452   - logL = ngam + alphap1 * nbgd - non - noff -
2453   - non * (std::log(nonpred) - std::log(non)) -
2454   - noff * (std::log(nbgd) - std::log(noff));
2455   -
2456   - // Compute derivatives
2457   - double f0 = alphat2 * noff - C;
2458   - double dCds = -alphap1;
2459   - double dDds = (C * dCds + 2.0 * alphap1 * alpha * noff) / D;
2460   - double dbds = (f0 / D - 1.0) / alphat2;
2461   - double d2bds2 = (-dCds / D - f0 / (D*D) * dDds) / alphat2;
2462   - double f1 = alphap1 - alpha*non/nonpred - noff/nbgd;
2463   - double f2 = nonpred * nonpred;
2464   - double dpds = 1.0 + alpha * dbds;
2465   - double f3 = non / f2 * dpds;
2466   - dlogLdsky = 1.0 - non / nonpred + dbds * f1;
2467   - d2logLdsky2 = f3 + d2bds2 * f1 +
2468   - dbds * (alpha * f3 + noff / (nbgd*nbgd) * dbds);
2469   -
2470   - } // endelse: general case
2471   -
2472   - // Return
2473   - return logL;
2474   -}
2475   -
2476   -
2477   -
2478   -/***********************************************************************//**
2479   - * @brief Evaluate log-likelihood function for On/Off analysis in the
2480   - * case of Poisson signal with measured Poisson background
2481   - *
2482   - * @param[in] models Models.
2483   - * @param[in,out] gradient Pointer to gradients.
2484   - * @param[in,out] curvature Pointer to curvature matrix.
2485   - * @param[in,out] npred Pointer to Npred value.
2486   - * @return Log-likelihood value.
2487   - *
2488   - * @exception GException::invalid_value
2489   - * There are no model parameters.
2490   - *
2491   - * Computes the log-likelihood value for the On/Off observation. The method
2492   - * loops over the energy bins to update the function value, its derivatives
2493   - * and the curvature matrix.
2494   - *
2495   - * The number of On counts \f$n_{\rm on}\f$ and Off counts \f$n_{\rm off}\f$
2496   - * are taken from the On and Off spectra, the expected number of gamma-ray
2497   - * events \f$\mu_s\f$ is computed from the spectral models of the relevant
2498   - * components in the model container (spatial and temporal components
2499   - * are ignored so far). See the N_gamma() method for details about the
2500   - * model computations.
2501   - *
2502   - * The estimated number of background counts \f$\mu_b\f$ is derived based
2503   - * on the measurement in the Off region by analytical minimization of
2504   - * the Poisson likelihood, i.e., it is treated as a nuisance parameter.
2505   - * See the wstat_value function for the derivation.
2506   - ***********************************************************************/
2507   -double GCTAOnOffObservation::likelihood_wstat(const GModels& models,
2508   - GVector* gradient,
2509   - GMatrixSparse* curvature,
2510   - double* npred) const
2511   -{
2512   - // Timing measurement
2513   - #if defined(G_LIKELIHOOD_DEBUG)
2514   - #ifdef _OPENMP
2515   - double t_start = omp_get_wtime();
2516   - #else
2517   - clock_t t_start = clock();
2518   - #endif
2519   - #endif
2520   -
2521   - // Initialise statistics
2522   - #if defined(G_LIKELIHOOD_DEBUG)
2523   - int n_bins = m_on_spec.size();
2524   - int n_used = 0;
2525   - int n_small_model = 0;
2526   - int n_zero_data = 0;
2527   - double sum_data = 0.0;
2528   - double sum_model = 0.0;
2529   - double init_npred = *npred;
2530   - #endif
2531   -
2532 2577 // Initialise log-likelihood value
2533   - double value = 0.0;
2534   -
2535   - // Get number of model parameters in model container
2536   - int npars = models.npars();
2537   -
2538   - // Create model gradient vectors for sky parameters
2539   - GVector sky_grad(npars);
2540   -
2541   - // Allocate working array
2542   - GVector colvar(npars);
2543   -
2544   - // Check that there is at least one parameter
2545   - if (npars > 0) {
2546   -
2547   - // Loop over all energy bins
2548   - for (int i = 0; i < m_on_spec.size(); ++i) {
2549   -
2550   - // Reinitialize working arrays
2551   - for (int j = 0; j < npars; ++j) {
2552   - sky_grad[j] = 0.0;
2553   - }
2554   -
2555   - // Get number of On and Off counts
2556   - double non = m_on_spec[i];
2557   - double noff = m_off_spec[i];
2558   -
2559   - // Get background scaling
2560   - double alpha = m_on_spec.backscal(i);
2561   -
2562   - // Get number of gamma and background events (and corresponding
2563   - // spectral model gradients)
2564   - double ngam = N_gamma(models, i, &sky_grad);
2565   -
2566   - // Initialise variables for likelihood computation
2567   - double nonpred = 0.0;
2568   - double nbgd = 0.0;
2569   - double dlogLdsky = 0.0;
2570   - double d2logLdsky2 = 0.0;
2571   -
2572   - // Perform likelihood profiling and derive likelihood value
2573   - double logL = wstat_value(non,noff,alpha,ngam,
2574   - nonpred,nbgd,dlogLdsky,d2logLdsky2);
2575   -
2576   - // Update global log-likelihood and Npred
2577   - value += logL;
2578   - *npred += nonpred;
2579   -
2580   - // Update statistics
2581   - #if defined(G_LIKELIHOOD_DEBUG)
2582   - n_used++;
2583   - sum_data += non;
2584   - sum_model += nonpred;
2585   - #endif
2586   -
2587   - // Loop over all parameters
2588   - for (int j = 0; j < npars; ++j) {
2589   -
2590   - // If spectral model for sky component is non-zero and
2591   - // non-infinite then handle sky component gradients and
2592   - // second derivatives including at least a sky component ...
2593   - if (sky_grad[j] != 0.0 && !gammalib::is_infinite(sky_grad[j])) {
2594   -
2595   - // Gradient
2596   - (*gradient)[j] += dlogLdsky * sky_grad[j];
2597   -
2598   - // Hessian (from first-order derivatives only)
2599   - for (int k = 0; k < npars; ++k) {
2600   -
2601   - // If spectral model for sky component is non-zero and
2602   - // non-infinite then we have the curvature element
2603   - // of a sky component
2604   - if (sky_grad[k] != 0.0 &&
2605   - !gammalib::is_infinite(sky_grad[k])) {
2606   - colvar[k] = sky_grad[j] * sky_grad[k] * d2logLdsky2;
2607   - }
2608   -
2609   - // ...... else if spectral model for background component
2610   - // or neither sky nor background
2611   - else {
2612   - colvar[k] = 0.0;
2613   - }
2614   -
2615   - } // endfor: Hessian computation
  2578 + double logL;
  2579 +
  2580 + // Precompute some values
  2581 + double alphap1 = alpha + 1.0;
  2582 + double alpharat = alpha / alphap1;
  2583 +
  2584 + // Calculate number of background events, profile likelihood value
  2585 + // and likelihood derivatives
  2586 + // Special case noff = 0
  2587 + if (noff == 0.0) {
  2588 +
  2589 + // Case A: non = 0. In this case nbgd < 0 hence we set nbgd = 0
  2590 + if (non == 0.0) {
  2591 + nbgd = 0;
  2592 + nonpred = ngam;
  2593 + logL = ngam;
  2594 + dlogLdsky = 1.0;
  2595 + d2logLdsky2 = 0.0;
  2596 + }
2616 2597  
2617   - // Update matrix
2618   - curvature->add_to_column(j, colvar);
  2598 + // Case B: nbgd is positive
  2599 + else if (ngam < non * alpharat) {
  2600 + nbgd = non / alphap1 - ngam / alpha;
  2601 + nonpred = ngam + alpha * nbgd;
  2602 + logL = -ngam / alpha - non * std::log(alpharat);
  2603 + dlogLdsky = -1.0/alpha;
  2604 + d2logLdsky2 = 0.0;
  2605 + }
2619 2606  
2620   - } // endif: spectral model is non-zero and non-infinite
  2607 + // Case C: nbgd is zero or negative, hence set nbgd = 0
  2608 + else {
  2609 + nbgd = 0;
  2610 + nonpred = ngam;
  2611 + logL = ngam + non * (std::log(non) - std::log(ngam) - 1.0);
  2612 + dlogLdsky = 1.0 - non / ngam;
  2613 + d2logLdsky2 = non / (ngam * ngam);
  2614 + }
2621 2615  
2622   - } // endfor: looped over all parameters for derivatives computation
  2616 + } // endif: noff = 0
2623 2617  
2624   - } // endfor: looped over energy bins
  2618 + // Special case non = 0
  2619 + else if (non == 0.0) {
  2620 + nbgd = 0;
  2621 + nonpred = ngam + alpharat * noff;
  2622 + logL = ngam + noff * std::log(alphap1);
  2623 + dlogLdsky = 1.0;
  2624 + d2logLdsky2 = 0.0;
  2625 + } // endif: non = 0
2625 2626  
2626   - } // endif: number of parameters was positive
2627   -
2628   - // ... else there are no parameters, so throw an exception
  2627 + // General case
2629 2628 else {
2630   - std::string msg ="No model parameter for the computation of the "
2631   - "likelihood in observation \""+this->name()+
2632   - "\" (ID "+this->id()+").\n";
2633   - throw GException::invalid_value(G_LIKELIHOOD_WSTAT, msg);
2634   - }
2635   -
2636   - // Optionally dump gradient and curvature matrix
2637   - #if defined(G_LIKELIHOOD_DEBUG)
2638   - std::cout << *gradient << std::endl;
2639   - std::cout << *curvature << std::endl;
2640   - #endif
2641 2629  
2642   - // Timing measurement
2643   - #if defined(G_LIKELIHOOD_DEBUG)
2644   - #ifdef _OPENMP
2645   - double t_elapse = omp_get_wtime()-t_start;
2646   - #else
2647   - double t_elapse = (double)(clock() - t_start) / (double)CLOCKS_PER_SEC;
2648   - #endif
2649   - std::cout << "GCTAOnOffObservation::optimizer::likelihood_wstat: CPU usage = "
2650   - << t_elapse << " sec" << std::endl;
2651   - #endif
  2630 + // Compute log-likelihood value
  2631 + double alphat2 = 2.0 * alpha;
  2632 + double C = alpha * (non + noff) - alphap1 * ngam;
  2633 + double D = std::sqrt(C*C + 4.0 * alpha * alphap1 * noff * ngam);
  2634 + nbgd = (C + D) / (alphat2 * alphap1);
  2635 + nonpred = ngam + alpha * nbgd;
  2636 + logL = ngam + alphap1 * nbgd - non - noff -
  2637 + non * (std::log(nonpred) - std::log(non)) -
  2638 + noff * (std::log(nbgd) - std::log(noff));
  2639 +
  2640 + // Compute derivatives
  2641 + double f0 = alphat2 * noff - C;
  2642 + double dCds = -alphap1;
  2643 + double dDds = (C * dCds + 2.0 * alphap1 * alpha * noff) / D;
  2644 + double dbds = (f0 / D - 1.0) / alphat2;
  2645 + double d2bds2 = (-dCds / D - f0 / (D*D) * dDds) / alphat2;
  2646 + double f1 = alphap1 - alpha*non/nonpred - noff/nbgd;
  2647 + double f2 = nonpred * nonpred;
  2648 + double dpds = 1.0 + alpha * dbds;
  2649 + double f3 = non / f2 * dpds;
  2650 + dlogLdsky = 1.0 - non / nonpred + dbds * f1;
  2651 + d2logLdsky2 = f3 + d2bds2 * f1 +
  2652 + dbds * (alpha * f3 + noff / (nbgd*nbgd) * dbds);
  2653 +
  2654 + } // endelse: general case
2652 2655  
2653 2656 // Return
2654   - return value;
  2657 + return logL;
2655 2658 }
2656 2659  
2657 2660  
... ...
pyext/GNdarray.i
  View file @e7028db
1 1 /***************************************************************************
2 2 * GNdarray.i - N-dimensional array class *
3 3 * ----------------------------------------------------------------------- *
4   - * copyright (C) 2016-2017 by Juergen Knoedlseder *
  4 + * copyright (C) 2016-2018 by Juergen Knoedlseder *
5 5 * ----------------------------------------------------------------------- *
6 6 * *
7 7 * This program is free software: you can redistribute it and/or modify *
... ...
pyext/GPha.i
  View file @e7028db
1 1 /***************************************************************************
2 2 * GPha.i - XSPEC Pulse Height Analyzer class *
3 3 * ----------------------------------------------------------------------- *
4   - * copyright (C) 2013-2017 by Juergen Knoedlseder *
  4 + * copyright (C) 2013-2018 by Juergen Knoedlseder *
5 5 * ----------------------------------------------------------------------- *
6 6 * *
7 7 * This program is free software: you can redistribute it and/or modify *
... ...
pyext/GSkyMap.i
  View file @e7028db
1 1 /***************************************************************************
2 2 * GSkyMap.i - Sky map class *
3 3 * ----------------------------------------------------------------------- *
4   - * copyright (C) 2010-2017 by Juergen Knoedlseder *
  4 + * copyright (C) 2010-2018 by Juergen Knoedlseder *
5 5 * ----------------------------------------------------------------------- *
6 6 * *
7 7 * This program is free software: you can redistribute it and/or modify *
... ... @@ -164,8 +164,10 @@ public:
164 164 int pix2inx(const GSkyPixel& pixel) const;
165 165 int dir2inx(const GSkyDir& dir) const;
166 166 GSkyPixel dir2pix(const GSkyDir& dir) const;
  167 + GNdarray counts(void);
167 168 double flux(const int& index, const int& map = 0) const;
168 169 double flux(const GSkyPixel& pixel, const int& map = 0) const;
  170 + GNdarray flux(void);
169 171 double solidangle(const int& index) const;
170 172 double solidangle(const GSkyPixel& pixel) const;
171 173 bool contains(const GSkyDir& dir) const;
... ... @@ -177,8 +179,6 @@ public:
177 179 const double* pixels(void) const;
178 180 GSkyMap extract(const int& map, const int& nmaps = 1) const;
179 181 void stack_maps(void);
180   - GNdarray total_counts(void);
181   - GNdarray total_flux(void);
182 182 void load(const GFilename& filename);
183 183 void save(const GFilename& filename,
184 184 const bool& clobber = false) const;
... ... @@ -259,7 +259,7 @@ public:
259 259 return log10(*self);
260 260 }
261 261 GSkyMap clip(const double& thresh) {
262   - return clip(*self,thresh);
  262 + return clip(*self,thresh);
263 263 }
264 264 GSkyMap abs() {
265 265 return abs(*self);
... ...
src/numerics/GNdarray.cpp
  View file @e7028db
1 1 /***************************************************************************
2 2 * GNdarray.cpp - N-dimensional array class *
3 3 * ----------------------------------------------------------------------- *
4   - * copyright (C) 2016-2017 by Juergen Knoedlseder *
  4 + * copyright (C) 2016-2018 by Juergen Knoedlseder *
5 5 * ----------------------------------------------------------------------- *
6 6 * *
7 7 * This program is free software: you can redistribute it and/or modify *
... ... @@ -333,6 +333,7 @@ GNdarray&amp; GNdarray::operator+=(const GNdarray&amp; array)
333 333 return *this;
334 334 }
335 335  
  336 +
336 337 /***********************************************************************//**
337 338 * @brief Unary multiplication operator
338 339 *
... ... @@ -355,6 +356,7 @@ GNdarray&amp; GNdarray::operator*=(const GNdarray&amp; array)
355 356 return *this;
356 357 }
357 358  
  359 +
358 360 /***********************************************************************//**
359 361 * @brief Unary division operator
360 362 *
... ... @@ -370,22 +372,24 @@ GNdarray&amp; GNdarray::operator/=(const GNdarray&amp; array)
370 372  
371 373 // Divide elements
372 374 for (int i = 0; i < m_data.size(); ++i) {
373   - // Only proceed if second array is non zero
374   - if (array.m_data[i] != 0.0){
375   - m_data[i] /= array.m_data[i];
376   - }
377   - // Otherwise set ratio to 0
378   - else {
379   - m_data[i] = 0.0;
380   - }
381   - }
  375 +
  376 + // Only proceed if second array is non zero
  377 + if (array.m_data[i] != 0.0){
  378 + m_data[i] /= array.m_data[i];
  379 + }
  380 +
  381 + // ... otherwise set ratio to 0
  382 + else {
  383 + m_data[i] = 0.0;
  384 + }
  385 +
  386 + } // endfor: looped over array elements
382 387  
383 388 // Return array
384 389 return *this;
385 390 }
386 391  
387 392  
388   -
389 393 /***********************************************************************//**
390 394 * @brief Unary subtraction operator
391 395 *
... ... @@ -1221,6 +1225,7 @@ GNdarray log(const GNdarray&amp; array)
1221 1225 return result;
1222 1226 }
1223 1227  
  1228 +
1224 1229 /***********************************************************************//**
1225 1230 * @brief Computes sign of array elements
1226 1231 *
... ... @@ -1234,26 +1239,28 @@ GNdarray sign(const GNdarray&amp; array)
1234 1239  
1235 1240 // Evaluate each array element
1236 1241 for (int i = 0; i < array.m_data.size(); ++i) {
1237   -
1238   - // Get array content
1239   - double content = array.m_data[i];
1240 1242  
1241   - // Initialise sign
1242   - double sign;
  1243 + // Get array content
  1244 + double content = array.m_data[i];
  1245 +
  1246 + // Initialise sign
  1247 + double sign;
1243 1248  
1244   - // Handle the 3 cases
1245   - if (content < 0.0) {
1246   - sign = -1.0;
1247   - }
1248   - else if (content > 0.0) {
1249   - sign = +1.0;
1250   - }
1251   - else {
1252   - sign = 0.0;
1253   - }
  1249 + // Handle the 3 cases
  1250 + if (content < 0.0) {
  1251 + sign = -1.0;
  1252 + }
  1253 + else if (content > 0.0) {
  1254 + sign = +1.0;
  1255 + }
  1256 + else {
  1257 + sign = 0.0;
  1258 + }
1254 1259  
  1260 + // Store sign
1255 1261 result.m_data[i] = sign;
1256   - }
  1262 +
  1263 + } // endfor: looped over all array elements
1257 1264  
1258 1265 // Return array
1259 1266 return result;
... ...
src/sky/GSkyMap.cpp
  View file @e7028db
1 1 /***************************************************************************
2 2 * GSkyMap.cpp - Sky map class *
3 3 * ----------------------------------------------------------------------- *
4   - * copyright (C) 2010-2017 by Juergen Knoedlseder *
  4 + * copyright (C) 2010-2018 by Juergen Knoedlseder *
5 5 * ----------------------------------------------------------------------- *
6 6 * *
7 7 * This program is free software: you can redistribute it and/or modify *
... ... @@ -1373,6 +1373,40 @@ GSkyPixel GSkyMap::dir2pix(const GSkyDir&amp; dir) const
1373 1373  
1374 1374  
1375 1375 /***********************************************************************//**
  1376 + * @brief Returns array with total number of counts for count maps
  1377 + *
  1378 + * @return Array of pixel counts sums.
  1379 + *
  1380 + * For a set of n count maps, the method returns a 1D array with n entries
  1381 + * that correspond to the sum of the pixel counts in each map.
  1382 + ***************************************************************************/
  1383 +GNdarray GSkyMap::counts(void) const
  1384 +{
  1385 + // Initialise output array
  1386 + GNdarray sum_array = GNdarray(m_num_maps);
  1387 +
  1388 + // Loop over maps and store sums in output array
  1389 + for (int i = 0; i < m_num_maps; ++i) {
  1390 +
  1391 + // Initialise sum for this map
  1392 + double sum = 0.0;
  1393 +
  1394 + // Compute map sum
  1395 + for (int k = 0; k < m_num_pixels; ++k) {
  1396 + sum += (*this)(k,i);
  1397 + }
  1398 +
  1399 + // Store map sum
  1400 + sum_array(i) = sum;
  1401 +
  1402 + } // endfor: looped over all maps
  1403 +
  1404 + // Return sum array
  1405 + return sum_array;
  1406 +}
  1407 +
  1408 +
  1409 +/***********************************************************************//**
1376 1410 * @brief Returns flux in pixel
1377 1411 *
1378 1412 * @param[in] index Pixel index [0,...,npix()-1].
... ... @@ -1581,6 +1615,41 @@ double GSkyMap::flux(const GSkyPixel&amp; pixel, const int&amp; map) const
1581 1615  
1582 1616  
1583 1617 /***********************************************************************//**
  1618 + * @brief Returns array with total flux for sky maps
  1619 + *
  1620 + * @return Array of pixel flux sums.
  1621 + *
  1622 + * For a set of n sky maps, the method returns a 1D array with n entries
  1623 + * that correspond to the sum of the pixel flux in each map.
  1624 + ***************************************************************************/
  1625 +GNdarray GSkyMap::flux(void) const
  1626 +{
  1627 + // Initialise output array
  1628 + GNdarray flux_array = GNdarray(m_num_maps);
  1629 +
  1630 + // Loop over maps and store sums in output array
  1631 + for (int i = 0; i < m_num_maps; ++i) {
  1632 +
  1633 + // Initialise flux for this map
  1634 + double flux = 0.0;
  1635 +
  1636 + // Compute map flux
  1637 + for (int k = 0; k < m_num_pixels; ++k) {
  1638 + double solidangle = this->solidangle(k);
  1639 + flux += (*this)(k,i) * solidangle;
  1640 + }
  1641 +
  1642 + // Store map flux
  1643 + flux_array(i) = flux;
  1644 +
  1645 + } // endfor: looped over all maps
  1646 +
  1647 + // Return flux array
  1648 + return flux_array;
  1649 +}
  1650 +
  1651 +
  1652 +/***********************************************************************//**
1584 1653 * @brief Returns solid angle of pixel
1585 1654 *
1586 1655 * @param[in] index Pixel index [0,...,npix()-1].
... ... @@ -1940,58 +2009,6 @@ void GSkyMap::stack_maps(void)
1940 2009 return;
1941 2010 }
1942 2011  
1943   -/***********************************************************************//**
1944   - * @brief Returns array with total number of counts for count maps
1945   - *
1946   - * @return Array of pixel counts sums.
1947   - *
1948   - * For a set of n count maps, the method returns a 1D array with n entries
1949   - * that correspond to the sum of the pixel counts in each map.
1950   - ***************************************************************************/
1951   -GNdarray GSkyMap::total_counts(void) const
1952   -{
1953   - // Initialise output array
1954   - GNdarray collapsed_array = GNdarray(m_num_maps);
1955   -
1956   - // Loop over maps and store sums in output array
1957   - for (int i = 0; i < m_num_maps; ++i) {
1958   - double sum = 0.0;
1959   - for (int k = 0; k < m_num_pixels; ++k) {
1960   - sum += (*this)(k,i);
1961   - }
1962   - collapsed_array(i) = sum;
1963   - }
1964   -
1965   - // Return pixel
1966   - return collapsed_array;
1967   -}
1968   -
1969   -/***********************************************************************//**
1970   - * @brief Returns array with total flux for sky maps
1971   - *
1972   - * @return Array of pixel flux sums.
1973   - *
1974   - * For a set of n sky maps, the method returns a 1D array with n entries
1975   - * that correspond to the sum of the pixel flux in each map.
1976   - ***************************************************************************/
1977   -GNdarray GSkyMap::total_flux(void) const
1978   -{
1979   - // Initialise output array
1980   - GNdarray collapsed_array = GNdarray(m_num_maps);
1981   -
1982   - // Loop over maps and store sums in output array
1983   - for (int i = 0; i < m_num_maps; ++i) {
1984   - double flux = 0.0;
1985   - for (int k = 0; k < m_num_pixels; ++k) {
1986   - double solidangle = this->solidangle(k);
1987   - flux += (*this)(k,i) * solidangle;
1988   - }
1989   - collapsed_array(i) = flux;
1990   - }
1991   -
1992   - // Return pixel
1993   - return collapsed_array;
1994   -}
1995 2012  
1996 2013 /***********************************************************************//**
1997 2014 * @brief Load skymap from FITS file.
... ... @@ -3296,44 +3313,6 @@ GSkyMap sqrt(const GSkyMap&amp; map)
3296 3313 return result;
3297 3314 }
3298 3315  
3299   -/***********************************************************************//**
3300   - * @brief Clips map at given value
3301   - *
3302   - * @param[in] map Sky map.
3303   - * @param[in] thresh Threshold value.
3304   - * @return Sky map containing the value in the original map if >= thresh,
3305   - * otherwise thresh
3306   - ***************************************************************************/
3307   -GSkyMap clip(const GSkyMap& map,const double& thresh)
3308   -{
3309   - // Initialise result vector
3310   - GSkyMap result(map);
3311   -
3312   - // Loop over all maps
3313   - for (int i = 0; i < map.nmaps(); ++i) {
3314   -
3315   - // Loop over all bins
3316   - for (int j = 0; j < map.npix(); ++j) {
3317   -
3318   - // Get the content from the bin
3319   - double content = map(j,i);
3320   -
3321   - // Ic content is below threshold, set to thresh
3322   - if (content < thresh) {
3323   - result(j,i) = thresh;
3324   - }
3325   -
3326   - else {
3327   - result(j,i) = content;
3328   - }
3329   -
3330   - } // endfor: looped over all bins
3331   -
3332   - } // endfor: looped over all maps
3333   -
3334   - // Return sky map
3335   - return result;
3336   -}
3337 3316  
3338 3317 /***********************************************************************//**
3339 3318 * @brief Computes the natural logarithm of sky map elements
... ... @@ -3489,3 +3468,44 @@ GSkyMap sign(const GSkyMap&amp; map)
3489 3468 // Return sky map
3490 3469 return result;
3491 3470 }
  3471 +
  3472 +
  3473 +/***********************************************************************//**
  3474 + * @brief Clips map at given value
  3475 + *
  3476 + * @param[in] map Sky map.
  3477 + * @param[in] thresh Threshold value.
  3478 + * @return Sky map containing the value in the original map if >= thresh,
  3479 + * otherwise thresh
  3480 + ***************************************************************************/
  3481 +GSkyMap clip(const GSkyMap& map, const double& thresh)
  3482 +{
  3483 + // Initialise result vector
  3484 + GSkyMap result(map);
  3485 +
  3486 + // Loop over all maps
  3487 + for (int i = 0; i < map.nmaps(); ++i) {
  3488 +
  3489 + // Loop over all bins
  3490 + for (int j = 0; j < map.npix(); ++j) {
  3491 +
  3492 + // Get the content from the bin
  3493 + double content = map(j,i);
  3494 +
  3495 + // If content is below threshold, set to thresh
  3496 + if (content < thresh) {
  3497 + result(j,i) = thresh;
  3498 + }
  3499 +
  3500 + // ... otherwise keep content
  3501 + else {
  3502 + result(j,i) = content;
  3503 + }
  3504 +
  3505 + } // endfor: looped over all bins
  3506 +
  3507 + } // endfor: looped over all maps
  3508 +
  3509 + // Return sky map
  3510 + return result;
  3511 +}
... ...
src/xspec/GPha.cpp
  View file @e7028db
1 1 /***************************************************************************
2 2 * GPha.cpp - XSPEC Pulse Height Analyzer class *
3 3 * ----------------------------------------------------------------------- *
4   - * copyright (C) 2013-2017 by Juergen Knoedlseder *
  4 + * copyright (C) 2013-2018 by Juergen Knoedlseder *
5 5 * ----------------------------------------------------------------------- *
6 6 * *
7 7 * This program is free software: you can redistribute it and/or modify *
... ... @@ -582,45 +582,51 @@ double GPha::counts(void) const
582 582 }
583 583  
584 584 /***********************************************************************//**
585   - * @brief Number of counts in spectrum as GNdArray
  585 + * @brief Get number of counts in spectrum as GNdarray
586 586 *
587   - * Returns the number of counts in the spectrum as a GNdArray.
  587 + * @return GNdarray of number of counts in spectrum.
  588 + *
  589 + * Returns the number of counts in the spectrum as GNdarray.
588 590 ***************************************************************************/
589 591 GNdarray GPha::counts_spectrum(void) const
590 592 {
591 593 // Initialise array
592   - int size = m_counts.size();
  594 + int size = m_counts.size();
593 595 GNdarray counts = GNdarray(size);
594   -
  596 +
595 597 // Compute content
596 598 for (int i = 0; i < size; ++i) {
597   - counts(i) = m_counts[i];
  599 + counts(i) = m_counts[i];
598 600 }
599 601  
600 602 // Return counts
601 603 return counts;
602 604 }
603 605  
  606 +
604 607 /***********************************************************************//**
605   - * @brief Number of counts in spectrum as GNdArray
  608 + * @brief Get background scaling factors as GNdarray
  609 + *
  610 + * @return GNdarray of background scaling factors.
606 611 *
607   - * Returns the number of counts in the spectrum as a GNdArray.
  612 + * Returns the background scaling factors as GNdarray.
608 613 ***************************************************************************/
609 614 GNdarray GPha::backscal_spectrum(void) const
610 615 {
611 616 // Initialise array
612   - int size = m_backscal.size();
  617 + int size = m_backscal.size();
613 618 GNdarray alpha = GNdarray(size);
614   -
  619 +
615 620 // Compute content
616 621 for (int i = 0; i < size; ++i) {
617   - alpha(i) = m_backscal[i];
  622 + alpha(i) = m_backscal[i];
618 623 }
619 624  
620   - // Return counts
  625 + // Return background scaling factors
621 626 return alpha;
622 627 }
623 628  
  629 +
624 630 /***********************************************************************//**
625 631 * @brief Fill spectrum with a value.
626 632 *
... ...
test/test_GSky.cpp
  View file @e7028db
1 1 /***************************************************************************
2 2 * test_GSky.cpp - Test sky module *
3 3 * ----------------------------------------------------------------------- *
4   - * copyright (C) 2010-2017 by Juergen Knoedlseder *
  4 + * copyright (C) 2010-2018 by Juergen Knoedlseder *
5 5 * ----------------------------------------------------------------------- *
6 6 * *
7 7 * This program is free software: you can redistribute it and/or modify *
... ... @@ -1056,17 +1056,17 @@ void TestGSky::test_GSkyMap(void)
1056 1056 test_value(map_stacked.nmaps(), 1, "Test stack_maps() method");
1057 1057  
1058 1058 // Test total counts computation
1059   - GNdarray counts_spectrum = map_src.total_counts();
  1059 + GNdarray counts_spectrum = map_src.counts();
1060 1060 double total_counts = 0.0;
1061 1061 for (int i = 0; i < counts_spectrum.size(); ++i) {
1062 1062 total_counts += counts_spectrum(i);
1063 1063 }
1064   - test_value(total_counts, total_src, 1.0e-3, "Test total_counts() method");
1065   - test_value(counts_spectrum.size(), map_src.nmaps(), "Test total_counts() method");
  1064 + test_value(total_counts, total_src, 1.0e-3, "Test counts() method");
  1065 + test_value(counts_spectrum.size(), map_src.nmaps(), "Test counts() method");
1066 1066  
1067 1067 // Test total flux computation
1068   - GNdarray flux_spectrum = map_src.total_flux();
1069   - test_value(flux_spectrum.size(), map_src.nmaps(), "Test total_flux() method");
  1068 + GNdarray flux_spectrum = map_src.flux();
  1069 + test_value(flux_spectrum.size(), map_src.nmaps(), "Test flux() method");
1070 1070  
1071 1071 // Test map number changing
1072 1072 GSkyMap map_more = map_src;
... ...