Johan Bregeon / gammalib

Commit 76da1e5bf983dd4a7e099eebc4fb78c7f9603897

Authored by Jürgen Knödlseder
1 parent e7b90fcd

Correct mc() method for spatial diffuse constant model. 

The formula for drawing random numbers on a sphere was not correct. This has now been corrected, and verify using a dedicated test script (test_model_spatial_const.py).
Inline Side-by-side
Showing 4 changed files with 73 additions and 6 deletions
ChangeLog
  View file @76da1e5
1   -2015-06-19 Juergen Knoedlseder <jurgen.knodlseder@irap.omp.eu>
  1 +2015-06-20 Juergen Knoedlseder <jurgen.knodlseder@irap.omp.eu>
2 2  
3 3 * Version 0.11.0 released
4 4 =========================
... ... @@ -6,6 +6,7 @@
6 6 Correct model gradient computation
7 7 Correct check of pixel containment in sky map (GSkymap::contains())
8 8 Do only search in users pfiles if PFILES environment variable is not set
  9 + Correct GModelSpatialDiffuseConst::mc() method (#1480)
9 10 Add is_valid() method to GCTARoi and write keywords only if ROI is valid (#1476)
10 11 Increased number of integration iterations for energy dispersion (#1472)
11 12 Add GOptimizerLM::logger() method, logger constructor takes pointer (#1471)
... ...
... ... @@ -2,7 +2,7 @@ New Features and Important Changes in GammaLib 0.11.0
2 2  
3 3 Author(s) : Juergen Knoedlseder
4 4  
5   -19 June 2015
  5 +20 June 2015
6 6  
7 7  
8 8 1. Introduction
... ... @@ -55,9 +55,9 @@ GModelSpatialDiffuseMap and GModelSpatialDiffuseCube. This assures that the
55 55 simulated photon distribution follows the intensity distribution of the skymap
56 56 as obtained using the sky direction access operator.
57 57  
58   -Integration of flux() and eflux() is now done using a fixed number of Romberg
59   -iterations instead of using the integration precision. This will assure that
60   -the log-likelihood function is smooth.
  58 +The algorithm used to draw random numbers in GModelSpatialDiffuseConst::mc()
  59 +was wrong (sky directions were not distributed isotropically on the sky). This
  60 +has now been fixed (#1480).
61 61  
62 62  
63 63 9. Numerics module
... ...
src/model/GModelSpatialDiffuseConst.cpp
  View file @76da1e5
... ... @@ -285,7 +285,7 @@ GSkyDir GModelSpatialDiffuseConst::mc(const GEnergy&amp; energy,
285 285 {
286 286 // Simulate Right Ascension and Declination
287 287 double ra = gammalib::twopi * ran.uniform();
288   - double dec = std::acos(1.0 - 2.0 * ran.uniform());
  288 + double dec = std::asin(2.0 * ran.uniform() - 1.0);
289 289  
290 290 // Set sky direction
291 291 GSkyDir dir;
... ...
test/test_model_spatial_const.py 0 → 100755
  View file @76da1e5
  1 +#! /usr/bin/env python
  2 +# ==========================================================================
  3 +# This script tests the GModelSpatialDiffuseConst model.
  4 +#
  5 +# Copyright (C) 2015 Juergen Knoedlseder
  6 +#
  7 +# This program is free software: you can redistribute it and/or modify
  8 +# it under the terms of the GNU General Public License as published by
  9 +# the Free Software Foundation, either version 3 of the License, or
  10 +# (at your option) any later version.
  11 +#
  12 +# This program is distributed in the hope that it will be useful,
  13 +# but WITHOUT ANY WARRANTY; without even the implied warranty of
  14 +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
  15 +# GNU General Public License for more details.
  16 +#
  17 +# You should have received a copy of the GNU General Public License
  18 +# along with this program. If not, see <http://www.gnu.org/licenses/>.
  19 +#
  20 +# ==========================================================================
  21 +import gammalib
  22 +#import math
  23 +
  24 +
  25 +# ============== #
  26 +# Evaluate model #
  27 +# ============== #
  28 +def test_mc(samples=100000):
  29 + """
  30 + This function tests the mc() method.
  31 + """
  32 + # Create model
  33 + model = gammalib.GModelSpatialDiffuseConst(1.0)
  34 +
  35 + # Create WCS map
  36 + map = gammalib.GSkymap("AIT", "GAL", 0.0, 0.0, -1.0, 1.0, 360, 180)
  37 +
  38 + # Set fixed parameters
  39 + energy = gammalib.GEnergy()
  40 + time = gammalib.GTime()
  41 +
  42 + # Allocate random number generator
  43 + ran = gammalib.GRan()
  44 +
  45 + # Fill map
  46 + for i in range(samples):
  47 + dir = model.mc(energy, time, ran)
  48 + pixel = map.dir2inx(dir)
  49 + map[pixel] += 1.0
  50 +
  51 + # Save map
  52 + map.save("test_model_spatial_const.fits", True)
  53 +
  54 + # Return
  55 + return
  56 +
  57 +
  58 +# ================ #
  59 +# Main entry point #
  60 +# ================ #
  61 +if __name__ == '__main__':
  62 + """
  63 + Main entry point
  64 + """
  65 + # Test Monte Carlo method
  66 + test_mc()
... ...