Alexander Ziegler / gammalib

Commit 5f79af820c16b65caaa596e60dd78f2a19d0a883

Authored by Jürgen Knödlseder
2 parents ef454174 f1a6d0d5

Merge branch '1725-GModelSpatialDiffuseMap-rejection' into integration 

# Conflicts:
#	ChangeLog
#	NEWS
#	README.md
Inline Side-by-side
Showing 8 changed files with 328 additions and 391 deletions
ChangeLog
  View file @5f79af8
... ... @@ -4,28 +4,29 @@
4 4 =========================
5 5  
6 6 Add GCTAPsfTable class (#1786)
7   - Use std::vector for naxes argument in GFitsImage (#1779).
8   - GCTACubePsf and GCTACubeEdisp now save 4D images.
9   - Change interpolation scheme in GCTABackground3D.
10   - Add GCsv::append() method (#1773).
11   - Add GTestSuite::test_value() string method (#1772).
12   - Add component setting methods to GModelSky (#1770).
13   - Add GApplication::logFileClose() method.
14   - Add GModel::read_attributes() and write_attributes() private methods (#1765).
15   - Also write test reports in case of an exception (#1746).
16   - Add GCTAObservation filename constructor (#1744).
17   - Add gammalib::number() function (#1743).
18   - Add GSkyMap::shape() and ndim() methods (#1733).
19   - Add energy dispersion cube for stacked analysis (#1721).
20   - Add GEnergy division operator.
  7 + Use std::vector for naxes argument in GFitsImage (#1779)
  8 + GCTACubePsf and GCTACubeEdisp now save 4D images
  9 + Change interpolation scheme in GCTABackground3D
  10 + Add GCsv::append() method (#1773)
  11 + Add GTestSuite::test_value() string method (#1772)
  12 + Add component setting methods to GModelSky (#1770)
  13 + Add GApplication::logFileClose() method
  14 + Add GModel::read_attributes() and write_attributes() private methods (#1765)
  15 + Also write test reports in case of an exception (#1746)
  16 + Add GCTAObservation filename constructor (#1744)
  17 + Add gammalib::number() function (#1743)
  18 + Add GSkyMap::shape() and ndim() methods (#1733)
  19 + Add GEnergy division operator
  20 + Use rejection method in GModelSpatialDiffuseMap::mc() (#1725)
21 21 Add GObservations::nobserved() method (#1723)
22 22 Use rejection method in GCTABackground3D::mc() (#1722)
23 23 Use rejection method in GModelSpatialDiffuseCube::mc() (#1722)
24   - Make GFilename::is_fits thread safe (#1717).
  24 + Add energy dispersion cube for stacked analysis (#1721)
  25 + Make GFilename::is_fits thread safe (#1717)
25 26 Revise GTime operators (#1714)
26 27 Add GSkyMap::publish() method to publish map on VO Hub (#1710)
27 28 Add GVOClient::publish() method to publish a FITS HDU (#1710)
28   - Remove GCTAOnOffObservations class (#1669).
  29 + Remove GCTAOnOffObservations class (#1669)
29 30 Handle FITS extensions in GSkyMap load() and save() methods (#1695)
30 31 Restrict Monte Carlo simulations in GModelSpatialDiffuseConst to simulation cone (#1674)
31 32 Use default GTime in GModelSpectral eval() and eval_gradients() methods (#1672)
... ... @@ -46,7 +47,7 @@
46 47 Use GFilename in GFits to store filename (#1649)
47 48 Load CTA GTI extension from the HDU name provided by the DSS keyword (#1598)
48 49 Support Fermi LAT Pass 8 response functions (#1517)
49   - Finalize implementation of CTA On/Off fitting (#1044).
  50 + Finalize implementation of CTA On/Off fitting (#1044)
50 51  
51 52  
52 53 2016-01-27 Juergen Knoedlseder <jurgen.knodlseder@irap.omp.eu>
... ...
... ... @@ -183,6 +183,9 @@ simulating sky directions, which makes the MC method identical to the model
183 183 evaluation method, using in both cases a log-log interpolation. The simulation
184 184 cone is now used for Monte Carlo simulations (#1722).
185 185  
  186 +Also the GModelSpatialDiffuseMap::mc() method now uses a rejection method for
  187 +simulating sky directions (#1725).
  188 +
186 189 The methods spatial(), spectral() and temporal() have been added to GModelSky
187 190 to allow setting of model components.
188 191  
... ...
include/GModelSpatialDiffuseMap.hpp
  View file @5f79af8
... ... @@ -110,11 +110,11 @@ protected:
110 110 bool m_has_normalize; //!< XML has normalize attribute
111 111  
112 112 // MC simulation cache
113   - mutable GSkyDir m_mc_centre; //!< Centre of MC cone
114   - mutable double m_mc_radius; //!< Radius of MC cone
115   - mutable double m_mc_norm; //!< Map normalization
116   - mutable std::vector<double> m_mc_cache; //!< Monte Carlo cache
117   - mutable std::vector<double> m_mc_max; //!< Monte Carlo maximum
  113 + mutable GSkyDir m_mc_centre; //!< Centre of MC cone
  114 + mutable double m_mc_radius; //!< Radius of MC cone
  115 + mutable double m_mc_one_minus_cosrad; //!< 1-cosine of radius
  116 + mutable double m_mc_norm; //!< Map normalization
  117 + mutable double m_mc_max; //!< Map maximum for MC
118 118 };
119 119  
120 120  
... ...
inst/cta/src/GCTACubeEdisp.cpp
  View file @5f79af8
... ... @@ -334,7 +334,7 @@ GCTACubeEdisp&amp; GCTACubeEdisp::operator=(const GCTACubeEdisp&amp; cube)
334 334 *
335 335 * @param[in] dir Coordinate of the true photon position.
336 336 * @param[in] migra Energy migration (reconstructed over true energy).
337   - * @param[in] energy True photon energy.
  337 + * @param[in] srcEng True photon energy.
338 338 * @return Energy dispersion (in units of MeV^-1)
339 339 *
340 340 * Returns the energy dispersion for a given energy migration (reconstructed
... ...
pyext/GModelSpatialDiffuseMap.i
  View file @5f79af8
... ... @@ -74,7 +74,7 @@ public:
74 74 void map(const GSkyMap& map);
75 75 bool normalize(void) const;
76 76 void set_mc_cone(const GSkyDir& centre,
77   - const double& radius) const;
  77 + const double& radius) const;
78 78 };
79 79  
80 80  
... ...
src/model/GModelSpatialDiffuseCube.cpp
  View file @5f79af8
... ... @@ -348,8 +348,9 @@ double GModelSpatialDiffuseCube::eval_gradients(const GPhoton&amp; photon) const
348 348 * @return Sky direction.
349 349 *
350 350 * @exception GException::invalid_value
351   - * No energy boundaries specified, or energy boundaries do not
352   - * cover the specified @p energy.
  351 + * No map cube defined.
  352 + * No energy boundaries defined.
  353 + * Simulation cone not defined or does not overlap with map cube.
353 354 *
354 355 * Returns a random sky direction according to the intensity distribution of
355 356 * the model sky map and the specified energy. The method uses a rejection
... ... @@ -359,102 +360,116 @@ GSkyDir GModelSpatialDiffuseCube::mc(const GEnergy&amp; energy,
359 360 const GTime& time,
360 361 GRan& ran) const
361 362 {
362   - // Allocate sky direction
363   - GSkyDir dir;
364   -
365 363 // Fetch cube
366 364 fetch_cube();
367 365  
368 366 // Determine number of skymap pixels
369 367 int npix = pixels();
370 368  
371   - // Continue only if there are skymap pixels
372   - if (npix > 0) {
373   -
374   - // If no energy boundaries are defined, throw an exception
375   - if (m_ebounds.size() < 1) {
376   - std::string msg = "The energy boundaries of the maps in the cube "
377   - "have not been defined. Maybe the map cube file "
378   - "is missing the \"ENERGIES\" extension which "
379   - "defines the energy of each map in the cube. "
380   - "Please provide the energy information.";
381   - throw GException::invalid_value(G_MC, msg);
382   - }
383   -
384   - // Set energy for interpolation
385   - m_logE.set_value(energy.log10MeV());
386   -
387   - // Compute maximum map value within simulation cone by log-log
388   - // interpolation of tha maximum values that are stored for each
389   - // map
390   - double max = 0.0;
391   - double max_left = m_mc_max[m_logE.inx_left()];
392   - double max_right = m_mc_max[m_logE.inx_right()];
393   - if (max_left > 0.0 && max_right > 0.0) {
394   - double max_log = m_logE.wgt_left() * std::log(max_left) +
395   - m_logE.wgt_right() * std::log(max_right);
396   - max = std::exp(max_log);
397   - }
398   - else if (max_left > 0.0) {
399   - max = max_left;
400   - }
401   - else if (max_right > 0.0) {
402   - max = max_right;
403   - }
404   -
405   - // Debug option: initialise counter
406   - #if defined(G_DEBUG_MC)
407   - int num_iterations = 0;
408   - #endif
  369 + // Throw an exception if there are no sky map pixels
  370 + if (npix <= 0) {
  371 + std::string msg = "No map cube defined. Please specify a valid map cube.";
  372 + throw GException::invalid_value(G_MC, msg);
  373 + }
409 374  
410   - // Get sky direction
411   - while (true) {
  375 + // Throw an exception if no energy boundaries are defined
  376 + if (m_ebounds.size() < 1) {
  377 + std::string msg = "The energy boundaries of the maps in the cube have "
  378 + "not been defined. Maybe the map cube file is missing "
  379 + "the \"ENERGIES\" extension which defines the energy "
  380 + "of each map in the cube. Please provide the energy "
  381 + "information.";
  382 + throw GException::invalid_value(G_MC, msg);
  383 + }
412 384  
413   - // Debug option: increment counter
414   - #if defined(G_DEBUG_MC)
415   - num_iterations++;
416   - #endif
  385 + // Set energy for interpolation
  386 + m_logE.set_value(energy.log10MeV());
  387 +
  388 + // Compute maximum map value within simulation cone by log-log interpolation
  389 + // of the maximum values that are stored for each map
  390 + double max = 0.0;
  391 + double max_left = m_mc_max[m_logE.inx_left()];
  392 + double max_right = m_mc_max[m_logE.inx_right()];
  393 + if (max_left > 0.0 && max_right > 0.0) {
  394 + double max_log = m_logE.wgt_left() * std::log(max_left) +
  395 + m_logE.wgt_right() * std::log(max_right);
  396 + max = std::exp(max_log);
  397 + }
  398 + else if (max_left > 0.0) {
  399 + max = max_left;
  400 + }
  401 + else if (max_right > 0.0) {
  402 + max = max_right;
  403 + }
417 404  
418   - // Simulate random sky direction within Monte Carlo simulation
419   - // cone
420   - double theta = std::acos(1.0 - ran.uniform() * m_mc_one_minus_cosrad) *
421   - gammalib::rad2deg;
422   - double phi = 360.0 * ran.uniform();
423   - dir = m_mc_centre;
424   - dir.rotate_deg(phi, theta);
425   -
426   - // Get map value at simulated sky direction
427   - double value = 0.0;
428   - double value_left = m_cube(dir, m_logE.inx_left());
429   - double value_right = m_cube(dir, m_logE.inx_right());
430   - if (value_left > 0.0 && value_right > 0.0) {
431   - double value_log = m_logE.wgt_left() * std::log(value_left) +
432   - m_logE.wgt_right() * std::log(value_right);
433   - value = std::exp(value_log);
434   - }
435   - else if (value_left > 0.0) {
436   - value = value_left;
437   - }
438   - else if (value_right > 0.0) {
439   - value = value_right;
440   - }
  405 + // Throw an exception if the maximum MC intensity is not positive. This
  406 + // can be the case because the simulation cone has not been defined or
  407 + // because it does not overlap with the sky map
  408 + if (max <= 0.0) {
  409 + std::string msg = "Simulation cone has not been defined or does not "
  410 + "overlap with the sky map. Please specify a valid "
  411 + "simulation cone.";
  412 + throw GException::invalid_value(G_MC, msg);
  413 + }
441 414  
442   - // Get uniform random number
443   - double uniform = ran.uniform() * max;
  415 + // Allocate sky direction
  416 + GSkyDir dir;
444 417  
445   - // Exit loop if we're not larger than the map value
446   - if (uniform <= value) {
447   - break;
448   - }
  418 + // Debug option: initialise counter
  419 + #if defined(G_DEBUG_MC)
  420 + int num_iterations = 0;
  421 + #endif
449 422  
450   - } // endwhile: loop until sky direction was accepted
  423 + // Get sky direction
  424 + while (true) {
451 425  
452   - // Debug option: log counter
  426 + // Debug option: increment counter
453 427 #if defined(G_DEBUG_MC)
454   - std::cout << num_iterations << " ";
  428 + num_iterations++;
455 429 #endif
456 430  
457   - } // endif: there were pixels in sky map
  431 + // Simulate random sky direction within Monte Carlo simulation cone
  432 + double theta = std::acos(1.0 - ran.uniform() * m_mc_one_minus_cosrad) *
  433 + gammalib::rad2deg;
  434 + double phi = 360.0 * ran.uniform();
  435 + dir = m_mc_centre;
  436 + dir.rotate_deg(phi, theta);
  437 +
  438 + // Get map value at simulated sky direction
  439 + double value = 0.0;
  440 + double value_left = m_cube(dir, m_logE.inx_left());
  441 + double value_right = m_cube(dir, m_logE.inx_right());
  442 + if (value_left > 0.0 && value_right > 0.0) {
  443 + double value_log = m_logE.wgt_left() * std::log(value_left) +
  444 + m_logE.wgt_right() * std::log(value_right);
  445 + value = std::exp(value_log);
  446 + }
  447 + else if (value_left > 0.0) {
  448 + value = value_left;
  449 + }
  450 + else if (value_right > 0.0) {
  451 + value = value_right;
  452 + }
  453 +
  454 + // If map value is non-positive then simulate a new sky direction
  455 + if (value <= 0.0) {
  456 + continue;
  457 + }
  458 +
  459 + // Get uniform random number
  460 + double uniform = ran.uniform() * max;
  461 +
  462 + // Exit loop if the random number is not larger than the map cube value
  463 + if (uniform <= value) {
  464 + break;
  465 + }
  466 +
  467 + } // endwhile: loop until sky direction was accepted
  468 +
  469 + // Debug option: log counter
  470 + #if defined(G_DEBUG_MC)
  471 + std::cout << num_iterations << " ";
  472 + #endif
458 473  
459 474 // Return sky direction
460 475 return dir;
... ... @@ -807,7 +822,7 @@ void GModelSpatialDiffuseCube::set_mc_cone(const GSkyDir&amp; centre,
807 822  
808 823 // Log maximum intensity and total flux for debugging
809 824 #if defined(G_DEBUG_MC_CACHE)
810   - std::cout << "GModelSpatialDiffuseCube::set_mc_cone:";
  825 + std::cout << "GModelSpatialDiffuseCube::set_mc_cone:" << std::endl;
811 826 std::cout << " Maximum map intensity:" << std::endl;
812 827 for (int i = 0; i < m_mc_max.size(); ++i) {
813 828 GEnergy energy;
... ...
src/model/GModelSpatialDiffuseMap.cpp
  View file @5f79af8
... ... @@ -42,6 +42,7 @@ const GModelSpatialDiffuseMap g_spatial_map_seed;
42 42 const GModelSpatialRegistry g_spatial_map_registry(&g_spatial_map_seed);
43 43  
44 44 /* __ Method name definitions ____________________________________________ */
  45 +#define G_MC "GModelSpatialDiffuseMap::mc(GEnergy&, GTime&, GRan&)"
45 46 #define G_READ "GModelSpatialDiffuseMap::read(GXmlElement&)"
46 47 #define G_WRITE "GModelSpatialDiffuseMap::write(GXmlElement&)"
47 48  
... ... @@ -50,8 +51,8 @@ const GModelSpatialRegistry g_spatial_map_registry(&amp;g_spatial_map_seed);
50 51 /* __ Coding definitions _________________________________________________ */
51 52  
52 53 /* __ Debug definitions __________________________________________________ */
53   -//#define G_DEBUG_CACHE //!< Debug Monte Carlo Cache computation
54   -//#define G_DEBUG_PREPARE //!< Debug map preparation
  54 +//#define G_DEBUG_MC //!< Debug MC method
  55 +//#define G_DEBUG_MC_CACHE //!< Debug MC cache
55 56  
56 57  
57 58 /*==========================================================================
... ... @@ -62,6 +63,8 @@ const GModelSpatialRegistry g_spatial_map_registry(&amp;g_spatial_map_seed);
62 63  
63 64 /***********************************************************************//**
64 65 * @brief Void constructor
  66 + *
  67 + * Constructs empty spatial map model.
65 68 ***************************************************************************/
66 69 GModelSpatialDiffuseMap::GModelSpatialDiffuseMap(void) :
67 70 GModelSpatialDiffuse()
... ... @@ -80,7 +83,7 @@ GModelSpatialDiffuseMap::GModelSpatialDiffuseMap(void) :
80 83 * @param[in] xml XML element.
81 84 *
82 85 * Constructs spatial map model by extracting information from an XML
83   - * element. See the method read() for more information about the expected
  86 + * element. See the read() method for more information about the expected
84 87 * structure of the XML element.
85 88 ***************************************************************************/
86 89 GModelSpatialDiffuseMap::GModelSpatialDiffuseMap(const GXmlElement& xml) :
... ... @@ -98,14 +101,15 @@ GModelSpatialDiffuseMap::GModelSpatialDiffuseMap(const GXmlElement&amp; xml) :
98 101  
99 102  
100 103 /***********************************************************************//**
101   - * @brief Filename constructor
  104 + * @brief File name constructor
102 105 *
103 106 * @param[in] filename File name.
104   - * @param[in] value Normalization factor (defaults to 1).
105   - * @param[in] normalize Normalize map (defaults to true).
  107 + * @param[in] value Normalization factor.
  108 + * @param[in] normalize Normalize map?
106 109 *
107   - * Constructs spatial map model by loading a skymap from the file specified
108   - * by @p filename and by setting the normalization @p value.
  110 + * Constructs spatial map model by loading a sky map from the file specified
  111 + * by @p filename and by setting the @p value by which the map will be
  112 + * multiplied (or normalized).
109 113 ***************************************************************************/
110 114 GModelSpatialDiffuseMap::GModelSpatialDiffuseMap(const GFilename& filename,
111 115 const double& value,
... ... @@ -115,13 +119,13 @@ GModelSpatialDiffuseMap::GModelSpatialDiffuseMap(const GFilename&amp; filename,
115 119 // Initialise members
116 120 init_members();
117 121  
118   - // Set normalization parameter
  122 + // Set normalisation parameter
119 123 m_value.value(value);
120 124  
121 125 // Set normalization flag
122 126 m_normalize = normalize;
123 127  
124   - // Load skymap
  128 + // Load sky map
125 129 load(filename);
126 130  
127 131 // Return
... ... @@ -133,11 +137,11 @@ GModelSpatialDiffuseMap::GModelSpatialDiffuseMap(const GFilename&amp; filename,
133 137 * @brief Sky map constructor
134 138 *
135 139 * @param[in] map Sky map.
136   - * @param[in] value Normalization factor (defaults to 1).
137   - * @param[in] normalize Normalize map (defaults to true).
  140 + * @param[in] value Normalization factor.
  141 + * @param[in] normalize Normalize map.
138 142 *
139 143 * Constructs spatial map model by setting the sky @p map and by setting the
140   - * normalization @p value.
  144 + * @p value by which the map will be multiplied (or normalized).
141 145 ***************************************************************************/
142 146 GModelSpatialDiffuseMap::GModelSpatialDiffuseMap(const GSkyMap& map,
143 147 const double& value,
... ... @@ -153,8 +157,10 @@ GModelSpatialDiffuseMap::GModelSpatialDiffuseMap(const GSkyMap&amp; map,
153 157 // Set normalization flag
154 158 m_normalize = normalize;
155 159  
156   - // Set and prepare skymap
  160 + // Set sky map
157 161 m_map = map;
  162 +
  163 + // Prepare sky map
158 164 prepare_map();
159 165  
160 166 // Return
... ... @@ -166,6 +172,8 @@ GModelSpatialDiffuseMap::GModelSpatialDiffuseMap(const GSkyMap&amp; map,
166 172 * @brief Copy constructor
167 173 *
168 174 * @param[in] model Spatial map model.
  175 + *
  176 + * Construct a spatial map model by copying another spatial map model.
169 177 ***************************************************************************/
170 178 GModelSpatialDiffuseMap::GModelSpatialDiffuseMap(const GModelSpatialDiffuseMap& model) :
171 179 GModelSpatialDiffuse(model)
... ... @@ -205,6 +213,8 @@ GModelSpatialDiffuseMap::~GModelSpatialDiffuseMap(void)
205 213 *
206 214 * @param[in] model Spatial map model.
207 215 * @return Spatial map model.
  216 + *
  217 + * Assigns a spatial map model to another spatial map model.
208 218 ***************************************************************************/
209 219 GModelSpatialDiffuseMap& GModelSpatialDiffuseMap::operator=(const GModelSpatialDiffuseMap& model)
210 220 {
... ... @@ -237,7 +247,7 @@ GModelSpatialDiffuseMap&amp; GModelSpatialDiffuseMap::operator=(const GModelSpatialD
237 247 ==========================================================================*/
238 248  
239 249 /***********************************************************************//**
240   - * @brief Clear diffuse map
  250 + * @brief Clear spatial map model
241 251 ***************************************************************************/
242 252 void GModelSpatialDiffuseMap::clear(void)
243 253 {
... ... @@ -257,9 +267,9 @@ void GModelSpatialDiffuseMap::clear(void)
257 267  
258 268  
259 269 /***********************************************************************//**
260   - * @brief Clone diffuse map
  270 + * @brief Clone spatial map model
261 271 *
262   - * @return Pointer to deep copy of diffuse map model.
  272 + * @return Pointer to deep copy of spatial map model.
263 273 ***************************************************************************/
264 274 GModelSpatialDiffuseMap* GModelSpatialDiffuseMap::clone(void) const
265 275 {
... ... @@ -269,14 +279,14 @@ GModelSpatialDiffuseMap* GModelSpatialDiffuseMap::clone(void) const
269 279  
270 280  
271 281 /***********************************************************************//**
272   - * @brief Return intensity of skymap
  282 + * @brief Return intensity of sky map
273 283 *
274 284 * @param[in] photon Incident photon.
275   - * @return Sky map intensity.
  285 + * @return Sky map intensity (\f$\mbox{ph cm}^{-2}\mbox{sr}^{-1}\mbox{s}^{-1}\f$)
276 286 *
277   - * Returns the intensity of the skymap at the specified sky direction
  287 + * Returns the intensity of the sky map at the specified sky direction
278 288 * multiplied by the normalization factor. If the sky direction falls outside
279   - * the skymap, an intensity of 0 is returned.
  289 + * the sky map, an intensity of 0 is returned.
280 290 ***************************************************************************/
281 291 double GModelSpatialDiffuseMap::eval(const GPhoton& photon) const
282 292 {
... ... @@ -289,15 +299,15 @@ double GModelSpatialDiffuseMap::eval(const GPhoton&amp; photon) const
289 299  
290 300  
291 301 /***********************************************************************//**
292   - * @brief Return intensity of skymap and gradient
  302 + * @brief Return intensity of sky map and gradient
293 303 *
294 304 * @param[in] photon Incident photon.
295   - * @return Sky map intensity.
  305 + * @return Sky map intensity (\f$\mbox{ph cm}^{-2}\mbox{sr}^{-1}\mbox{s}^{-1}\f$)
296 306 *
297   - * Returns the intensity of the skymap at the specified sky direction
  307 + * Returns the intensity of the sky map at the specified sky direction
298 308 * multiplied by the normalization factor. The method also sets the gradient
299 309 * with respect to the normalization factor. If the sky direction falls
300   - * outside the skymap, an intensity of 0 is returned.
  310 + * outside the sky map, an intensity of 0 is returned.
301 311 ***************************************************************************/
302 312 double GModelSpatialDiffuseMap::eval_gradients(const GPhoton& photon) const
303 313 {
... ... @@ -323,121 +333,81 @@ double GModelSpatialDiffuseMap::eval_gradients(const GPhoton&amp; photon) const
323 333 * @param[in,out] ran Random number generator.
324 334 * @return Sky direction.
325 335 *
326   - * Draws a random sky direction from the intensity distribution of the model
327   - * sky map.
328   - *
329   - * The method makes use of a cache array that contains the normalized
330   - * cumulative flux values of the skymap. Using a uniform random number, this
331   - * cache array is scanned using a bi-section method to determine the skymap
332   - * pixel for which the position should be returned.
333   - *
334   - * Within that pixel, a rejection method is used to draw a sky direction
335   - * that follows the intensity distribution that is obtained when calling the
336   - * interpolation operator. This assures that even for coarse binning of the
337   - * sky map the simulation corresponds to the model.
  336 + * @exception GException::invalid_value
  337 + * Simulation cone not defined or does not overlap with sky map.
  338 + * No sky map defined.
338 339 *
339   - * Note that the set_mc_cone() method needs to be called before this method
340   - * is called the first time.
  340 + * Returns a random sky direction according to the intensity distribution of
  341 + * the model sky map. The method uses a rejection method to determine the sky
  342 + * direction.
341 343 ***************************************************************************/
342 344 GSkyDir GModelSpatialDiffuseMap::mc(const GEnergy& energy,
343 345 const GTime& time,
344 346 GRan& ran) const
345 347 {
346   - // Allocate sky direction
347   - GSkyDir dir;
  348 + // Throw an exception if the maximum MC intensity is not positive. This
  349 + // can be the case because the simulation cone has not been defined or
  350 + // because it does not overlap with the sky map
  351 + if (m_mc_max <= 0.0) {
  352 + std::string msg = "Simulation cone has not been defined or does not "
  353 + "overlap with the sky map. Please specify a valid "
  354 + "simulation cone.";
  355 + throw GException::invalid_value(G_MC, msg);
  356 + }
348 357  
349 358 // Determine number of skymap pixels
350 359 int npix = m_map.npix();
351 360  
352   - // Continue only if there are skymap pixels
353   - if (npix > 0) {
354   -
355   - // Get pixel index from CDF
356   - int index = ran.cdf(m_mc_cache);
357   -
358   - // Convert sky map index to sky map pixel
359   - GSkyPixel pixel = m_map.inx2pix(index);
360   -
361   - // If we have a 2D pixel then randomize pixel values and convert them
362   - // into a sky direction
363   - if (pixel.is_2D()) {
364   -
365   - // Use a rejection method to find a sky direction
366   - while (true) {
367   -
368   - // Draw random sky pixel
369   - GSkyPixel test(pixel.x() + ran.uniform() - 0.5,
370   - pixel.y() + ran.uniform() - 0.5);
371   -
372   - // Derive sky direction
373   - dir = m_map.pix2dir(test);
374   -
375   - // Get map value at that sky direction
376   - double value = m_map(dir);
377   -
378   - // Get uniform random number up to the maximum
379   - double uniform = ran.uniform() * m_mc_max[index];
380   -
381   - // Exit loop if we're not larger than the map value
382   - if (uniform <= value) {
383   - break;
384   - }
385   -
386   - } // endwhile: rejection method
387   -
388   - } // endif: had a 2D pixel
389   -
390   - // ... otherwise convert pixel into sky direction and randomize
391   - // position. We use a while loop to randomize the position within a
392   - // circle that enclosed the pixel and retain only directions that
393   - // result in the same pixel index and that are compatible with the
394   - // density distribution.
395   - else {
396   -
397   - // Get pointer on HealPix projection
398   - const GHealpix* healpix = static_cast<const GHealpix*>(m_map.projection());
399   -
400   - // Get enclosing radius
401   - double radius = healpix->max_pixrad();
402   -
403   - // Initialize pixel centre
404   - dir = m_map.pix2dir(pixel);
405   -
406   - // Get randomized pixel
407   - GSkyDir randomized_dir;
408   - double cosrad = std::cos(radius);
409   - while (true) {
410   -
411   - // Get randomized sky direction
412   - randomized_dir = dir;
413   - double theta = std::acos(1.0 - ran.uniform() * (1.0 - cosrad)) * gammalib::rad2deg;
414   - double phi = 360.0 * ran.uniform();
415   - randomized_dir.rotate_deg(phi, theta);
416   -
417   - // Skip if we are not in the actual pixel
418   - if (m_map.dir2inx(randomized_dir) != index) {
419   - continue;
420   - }
421   -
422   - // Get map value at that sky direction
423   - double value = m_map(randomized_dir);
  361 + // Throw an exception if there are no sky map pixels
  362 + if (npix <= 0) {
  363 + std::string msg = "No sky map defined. Please specify a valid sky map.";
  364 + throw GException::invalid_value(G_MC, msg);
  365 + }
424 366  
425   - // Get uniform random number up to the maximum
426   - double uniform = ran.uniform() * m_mc_max[index];
  367 + // Allocate sky direction
  368 + GSkyDir dir;
427 369  
428   - // Exit loop if we're not larger than the map value
429   - if (uniform <= value) {
430   - break;
431   - }
  370 + // Debug option: initialise counter
  371 + #if defined(G_DEBUG_MC)
  372 + int num_iterations = 0;
  373 + #endif
  374 +
  375 + // Get sky direction
  376 + while (true) {
  377 +
  378 + // Debug option: increment counter
  379 + #if defined(G_DEBUG_MC)
  380 + num_iterations++;
  381 + #endif
  382 +
  383 + // Simulate random sky direction within Monte Carlo simulation cone
  384 + double theta = std::acos(1.0 - ran.uniform() * m_mc_one_minus_cosrad) *
  385 + gammalib::rad2deg;
  386 + double phi = 360.0 * ran.uniform();
  387 + dir = m_mc_centre;
  388 + dir.rotate_deg(phi, theta);
  389 +
  390 + // Get map value at simulated sky direction. If the map value is non-
  391 + // positive then simulate a new sky direction.
  392 + double value = m_map(dir);
  393 + if (value <= 0.0) {
  394 + continue;
  395 + }
432 396  
433   - } // endwhile
  397 + // Get uniform random number
  398 + double uniform = ran.uniform() * m_mc_max;
434 399  
435   - // Store randomize sky position
436   - dir = randomized_dir;
  400 + // Exit loop if the random number is not larger than the sky map value
  401 + if (uniform <= value) {
  402 + break;
  403 + }
437 404  
438   - } // endelse: we had a HealPix map
  405 + } // endwhile: loop until sky direction was accepted
439 406  
440   - } // endif: there were pixels in sky map
  407 + // Debug option: log counter
  408 + #if defined(G_DEBUG_MC)
  409 + std::cout << num_iterations << " ";
  410 + #endif
441 411  
442 412 // Return sky direction
443 413 return dir;
... ... @@ -517,7 +487,7 @@ bool GModelSpatialDiffuseMap::contains(const GSkyDir&amp; dir,
517 487 * Read the skymap information from an XML element. The XML element is
518 488 * required to have 1 parameter named either "Normalization" or "Prefactor".
519 489 *
520   - * If the attribute @a normalize="0" or @ normalize="false" is present the
  490 + * If the attribute @a normalize="0" or @a normalize="false" is present the
521 491 * diffuse map will not be normalised to unity flux upon loading.
522 492 ***************************************************************************/
523 493 void GModelSpatialDiffuseMap::read(const GXmlElement& xml)
... ... @@ -643,13 +613,9 @@ void GModelSpatialDiffuseMap::write(GXmlElement&amp; xml) const
643 613 * @param[in] radius Simulation cone radius (degrees).
644 614 *
645 615 * Sets the simulation cone centre and radius that defines the directions
646   - * that will be simulated using the mc() method.
647   - *
648   - * The method initialises a cache for Monte Carlo sampling of the skymap.
649   - * This Monte Carlo cache consists of a linear array that maps a value
650   - * between 0 and 1 into the skymap pixel. A second array contains the
651   - * maximum expected value for each pixel which is also used in Monte Carlo
652   - * sampling.
  616 + * that will be simulated using the mc() method and pre-computes the maximum
  617 + * intensity and the spatially integrated flux of the map within the
  618 + * simulation cone region.
653 619 ***************************************************************************/
654 620 void GModelSpatialDiffuseMap::set_mc_cone(const GSkyDir& centre,
655 621 const double& radius) const
... ... @@ -661,9 +627,11 @@ void GModelSpatialDiffuseMap::set_mc_cone(const GSkyDir&amp; centre,
661 627 m_mc_centre = centre;
662 628 m_mc_radius = radius;
663 629  
664   - // Initialise cache
665   - m_mc_cache.clear();
666   - m_mc_max.clear();
  630 + // Pre-compute 1 - cosine of radius
  631 + m_mc_one_minus_cosrad = 1.0 - std::cos(m_mc_radius*gammalib::deg2rad);
  632 +
  633 + // Initialise map maximum and normalisation
  634 + m_mc_max = 0.0;
667 635 m_mc_norm = 0.0;
668 636  
669 637 // Determine number of skymap pixels
... ... @@ -672,63 +640,36 @@ void GModelSpatialDiffuseMap::set_mc_cone(const GSkyDir&amp; centre,
672 640 // Continue only if there are pixels
673 641 if (npix > 0) {
674 642  
675   - // Reserve space for all pixels in cache
676   - m_mc_cache.reserve(npix+1);
677   - m_mc_max.reserve(npix);
678   -
679   - // Set first cache value to 0
680   - m_mc_cache.push_back(0.0);
681   -
682   - // Initialise cache with cumulative pixel fluxes and compute
683   - // total flux in skymap for normalization. Negative pixels are
684   - // excluded from the cumulative map. Invalid pixels are also
685   - // filtered.
  643 + // Compute flux and maximum map intensity within the simulation cone
686 644 double sum = 0.0;
687 645 double sum_map = 0.0;
688 646 for (int i = 0; i < npix; ++i) {
689 647  
690   - // Derive effective pixel radius from half opening angle
691   - // that corresponds to the pixel's solid angle. For security,
692   - // the radius is enhanced by 50%.
693   - double pixel_radius =
694   - std::acos(1.0 - m_map.solidangle(i)/gammalib::twopi) *
695   - gammalib::rad2deg * 1.5;
696   -
697   - // Add up flux with simulation cone radius + effective pixel
698   - // radius. The effective pixel radius is added to make sure
699   - // that all pixels that overlap with the simulation cone are
700   - // taken into account. There is no problem of having even
701   - // pixels outside the simulation cone taken into account as
702   - // long as the mc() method has an explicit test of whether a
703   - // simulated event is contained in the simulation cone.
704   - double flux = m_map.flux(i);
  648 + // Get map flux, intensity and distance from MC cone centre
  649 + double flux = m_map.flux(i);
  650 + double intensity = m_map(i);
  651 + double distance = centre.dist_deg(m_map.pix2dir(i));
  652 +
  653 + // Add flux if positive
705 654 if (flux > 0.0) {
706   - double distance = centre.dist_deg(m_map.pix2dir(i));
707   - if (distance <= radius+pixel_radius) {
708   - sum += flux;
  655 + if (distance <= radius) {
  656 + sum += flux; // flux within simulation cone
709 657 }
710   - sum_map += flux; // sum up total flux in map
  658 + sum_map += flux; // total flux
711 659 }
712   -
713   - // Push back flux
714   - m_mc_cache.push_back(sum); // units: ph/cm2/s/MeV
715   -
716   - } // endfor: looped over pixels
717   -
718   - // Normalize fluxes in the cache so that the values in the cache
719   - // run from 0 to 1.
720   - if (sum > 0.0) {
721   - for (int i = 0; i < npix; ++i) {
722   - m_mc_cache[i] /= sum;
  660 +
  661 + // Update maximum intensity
  662 + if (distance <= radius) {
  663 + if (intensity > m_mc_max) {
  664 + m_mc_max = intensity;
  665 + }
723 666 }
724   - }
725 667  
726   - // Make sure that last pixel in the cache is >1
727   - m_mc_cache[npix] = 1.0001;
  668 + } // endfor: looped over all map pixels
728 669  
729 670 // Set the normalization factor for the MC simulations. In case
730   - // that the map is normalized, this is the fraction of the flux
731   - // that is comprised within the simulation cone. For non-normalized
  671 + // that the map is normalised, this is the fraction of the flux
  672 + // that is comprised within the simulation cone. For non-normalised
732 673 // maps, this is simply the flux comprised within the simulation
733 674 // cone.
734 675 if (sum_map > 0.0) {
... ... @@ -740,79 +681,18 @@ void GModelSpatialDiffuseMap::set_mc_cone(const GSkyDir&amp; centre,
740 681 }
741 682 }
742 683  
743   - // Do we have a HealPix map?
744   - if (m_map.projection()->code() == "HPX") {
745   -
746   - // Get pointer on HealPix projection
747   - const GHealpix* healpix =
748   - static_cast<const GHealpix*>(m_map.projection());
749   -
750   - // Compute maximum value that may occur from bilinear
751   - // interpolation within this pixel and push this value on the
752   - // stack. We do this by checking values of all neighbours.
753   - for (int i = 0; i < npix; ++i) {
754   -
755   - // Get neighbours
756   - std::vector<int> neighbours = healpix->neighbours(i);
757   -
758   - // Loop over neighbours
759   - double max = m_map(i);
760   - for (int j = 0; j < neighbours.size(); ++j) {
761   - if (neighbours[j] != -1) {
762   - double value = m_map(neighbours[j]);
763   - if (value > max) {
764   - max = value;
765   - }
766   - }
767   - }
768   -
769   - // Store maximum
770   - m_mc_max.push_back(max);
771   -
772   - } // endfor: looped over pixels
773   -
774   - } // endif: Healpix projection
775   -
776   - // ... no, then we have a WCS map
777   - else {
778   -
779   - // Compute maximum value that may occur from bilinear
780   - // interpolation within this pixel and push this value on the
781   - // stack. We do this by checking the map values at the corners
782   - // and the centre of each edge.
783   - for (int i = 0; i < npix; ++i) {
784   - GSkyPixel pixel = m_map.inx2pix(i);
785   - double max = m_map(pixel);
786   - for (int ix = -1; ix < 2; ++ix) {
787   - for (int iy = -1; iy < 2; ++iy) {
788   - if (ix != 0 && iy != 0) {
789   - GSkyPixel edge(pixel.x()+ix*0.5, pixel.y()+iy*0.5);
790   - if (m_map.contains(edge)) {
791   - GSkyDir dir = m_map.pix2dir(edge);
792   - double value = m_map(dir);
793   - if (value > max) {
794   - max = value;
795   - }
796   - }
797   - }
798   - }
799   - }
800   - m_mc_max.push_back(max);
801   - }
802   -
803   - } // endelse: computed maximum pixel value
804   -
805   - // Dump cache values for debugging
806   - #if defined(G_DEBUG_CACHE)
807   - std::cout << "GModelSpatialDiffuseMap::set_mc_cone: cache";
808   - std::cout << std::endl;
809   - for (int i = 0; i < m_mc_cache.size(); ++i) {
810   - std::cout << "i=" << i;
811   - std::cout << " c=" << m_mc_cache[i] << std::endl;
812   - }
  684 + // Log maximum intensity and total flux for debugging
  685 + #if defined(G_DEBUG_MC_CACHE)
  686 + std::cout << "GModelSpatialDiffuseMap::set_mc_cone:" << std::endl;
  687 + std::cout << " Maximum map intensity:";
  688 + std::cout << m_mc_max << " ph/cm2/s/sr" << std::endl;
  689 + std::cout << " Spatially integrated flux:" << std::endl;
  690 + std::cout << sum << " ph/cm2/s" << std::endl;
  691 + std::cout << " Map normalisation:" << std::endl;
  692 + std::cout << m_mc_norm << " ph/cm2/s" << std::endl;
813 693 #endif
814 694  
815   - } // endif: there were cube pixels
  695 + } // endif: there were map pixels
816 696  
817 697 } // endif: simulation cone has changed
818 698  
... ... @@ -824,7 +704,7 @@ void GModelSpatialDiffuseMap::set_mc_cone(const GSkyDir&amp; centre,
824 704 /***********************************************************************//**
825 705 * @brief Print map information
826 706 *
827   - * @param[in] chatter Chattiness (defaults to NORMAL).
  707 + * @param[in] chatter Chattiness.
828 708 * @return String with diffuse map model information.
829 709 ***************************************************************************/
830 710 std::string GModelSpatialDiffuseMap::print(const GChatter& chatter) const
... ... @@ -882,7 +762,7 @@ void GModelSpatialDiffuseMap::load(const GFilename&amp; filename)
882 762 // variables
883 763 m_filename = filename;
884 764  
885   - // Load skymap
  765 + // Load sky map
886 766 m_map.load(m_filename);
887 767  
888 768 // Prepare sky map
... ... @@ -928,10 +808,10 @@ void GModelSpatialDiffuseMap::init_members(void)
928 808  
929 809 // Initialise MC cache
930 810 m_mc_centre.clear();
931   - m_mc_radius = -1.0; //!< Signal initialisation
932   - m_mc_norm = 0.0;
933   - m_mc_cache.clear();
934   - m_mc_max.clear();
  811 + m_mc_radius = -1.0; //!< Signal initialisation
  812 + m_mc_one_minus_cosrad = 1.0;
  813 + m_mc_norm = 0.0;
  814 + m_mc_max = 0.0;
935 815  
936 816 // Return
937 817 return;
... ... @@ -955,11 +835,11 @@ void GModelSpatialDiffuseMap::copy_members(const GModelSpatialDiffuseMap&amp; model)
955 835 m_radius = model.m_radius;
956 836  
957 837 // Copy MC cache
958   - m_mc_centre = model.m_mc_centre;
959   - m_mc_radius = model.m_mc_radius;
960   - m_mc_norm = model.m_mc_norm;
961   - m_mc_cache = model.m_mc_cache;
962   - m_mc_max = model.m_mc_max;
  838 + m_mc_centre = model.m_mc_centre;
  839 + m_mc_radius = model.m_mc_radius;
  840 + m_mc_one_minus_cosrad = model.m_mc_one_minus_cosrad;
  841 + m_mc_norm = model.m_mc_norm;
  842 + m_mc_max = model.m_mc_max;
963 843  
964 844 // Set parameter pointer(s)
965 845 m_pars.clear();
... ... @@ -983,12 +863,14 @@ void GModelSpatialDiffuseMap::free_members(void)
983 863 /***********************************************************************//**
984 864 * @brief Prepare sky map after loading
985 865 *
986   - * Prepares a sky map after loading. The map is normalized so that the total
987   - * flux in the map amounts to 1 ph/cm2/s. Negative skymap pixels are set to
988   - * zero intensity.
  866 + * Prepares a sky map after loading. Negative, infinite or undefined skymap
  867 + * pixels are set to zero intensity. The method also determine the centre
  868 + * and radius of a circle enclosing the map, and the Monte Carlo simulation
  869 + * cone is set to this circle.
989 870 *
990   - * Note that if the GSkyMap object contains multiple maps, only the first
991   - * map is used.
  871 + * If normalize() is true, the map is furthermore normalised so that the
  872 + * total flux in the map amounts to 1 ph/cm2/s. Negative skymap pixels are
  873 + * set to zero intensity.
992 874 ***************************************************************************/
993 875 void GModelSpatialDiffuseMap::prepare_map(void)
994 876 {
... ... @@ -1049,6 +931,9 @@ void GModelSpatialDiffuseMap::prepare_map(void)
1049 931  
1050 932 } // endelse: computed map centre and radius
1051 933  
  934 + // Set simulation cone
  935 + set_mc_cone(m_centre, m_radius);
  936 +
1052 937 } // endif: there were skymap pixels
1053 938  
1054 939 // Return
... ...
test/test_model_spatial_map.py
  View file @5f79af8
1 1 #! /usr/bin/env python
2 2 # ==========================================================================
3   -# This script tests the GModelSpatialMap spatial map model.
  3 +# This script tests the spatial map model GModelSpatialDiffuseMap.
4 4 #
5   -# Copyright (C) 2012 Juergen Knoedlseder
  5 +# Copyright (C) 2012-2016 Juergen Knoedlseder
6 6 #
7 7 # This program is free software: you can redistribute it and/or modify
8 8 # it under the terms of the GNU General Public License as published by
... ... @@ -18,29 +18,62 @@
18 18 # along with this program. If not, see <http://www.gnu.org/licenses/>.
19 19 #
20 20 # ==========================================================================
21   -from gammalib import *
22 21 import math
  22 +import gammalib
23 23  
24 24  
25 25 # ============== #
26 26 # Evaluate model #
27 27 # ============== #
28   -def eval_model(model, binsz=0.01, nbins=1000):
  28 +def eval_model(model, binsz=0.02, nbins=200, ntrials=1000):
29 29 """
30   - This function uses the eval method of the GModelSpatialMap class to
  30 + This function uses the eval method of the GModelSpatialDiffuseMap class to
31 31 interpolate the skymap into a finer grid.
32 32 """
  33 + # Define centre and radius
  34 + centre = gammalib.GSkyDir()
  35 + centre.radec_deg(201.3651, -43.0191)
  36 + radius = 1.5
  37 +
33 38 # Create CAR skymap
34   - map = GSkyMap("CAR", "CEL", 201.3651, -43.0191, binsz, binsz, nbins, nbins)
  39 + map = gammalib.GSkyMap("CAR", "CEL", 201.3651, -43.0191, binsz, binsz, nbins, nbins)
35 40  
36 41 # Fill map
  42 + sum = 0.0
  43 + for i in range(map.npix()):
  44 + dir = map.inx2dir(i)
  45 + if centre.dist_deg(dir) <= radius:
  46 + photon = gammalib.GPhoton(dir, gammalib.GEnergy(), gammalib.GTime())
  47 + intensity = model.eval(photon) * map.solidangle(i)
  48 + map[i] = intensity
  49 + sum += intensity
  50 +
  51 + # Normalize map to number of trials
  52 + norm = model.mc_norm(centre, radius)
  53 + print(norm)
  54 + print(sum)
  55 + print(norm/sum)
37 56 for i in range(map.npix()):
38   - dir = map.pix2dir(i)
39   - intensity = model.eval(dir)
40   - map[i] = intensity
  57 + map[i] = -map[i] * float(ntrials)/norm
  58 +
  59 + # Add MC simulated sky directions
  60 + ran = gammalib.GRan()
  61 + for i in range(ntrials):
  62 + dir = model.mc(gammalib.GEnergy(), gammalib.GTime(), ran)
  63 + inx = map.dir2inx(dir)
  64 + if inx >= 0 and inx < map.npix():
  65 + map[inx] += 1.0
  66 + else:
  67 + print('Pixel %d out of map boundary' % inx)
  68 +
  69 + # Determine map residuals
  70 + residual = 0.0
  71 + for i in range(map.npix()):
  72 + residual += map[i]
  73 + print('Residual: %e (%8.5f)' % (residual, residual/float(ntrials)))
41 74  
42 75 # Save map
43   - map.save("test_model_spatial_map.fits")
  76 + map.save('test_model_spatial_map.fits', True)
44 77  
45 78 # Return
46 79 return
... ... @@ -50,12 +83,12 @@ def eval_model(model, binsz=0.01, nbins=1000):
50 83 # Main entry point #
51 84 # ================ #
52 85 if __name__ == '__main__':
53   - """
54   - Main entry point
55   - """
  86 +
56 87 # Allocate spatial map model
57   - model = GModelSpatialMap("data/cena_lobes_parkes.fits")
58   - print model
  88 + model = gammalib.GModelSpatialDiffuseMap('data/cena_lobes_parkes.fits')
  89 + print(model)
59 90  
60 91 # Evaluate skymap
61   - eval_model(model)
  92 + #eval_model(model)
  93 + #eval_model(model, ntrials=10000000)
  94 + eval_model(model, ntrials=10000000)
... ...