Luigi Tibaldo / gammalib

Commit c1d193c1ea325a1c7a510a9f25d392ced0cdc435

Authored by Jürgen Knödlseder
1 parent 38f25250

Fix GModelSpatialDiffuseCube pickeling issue (#3179)

Inline Side-by-side
Showing 11 changed files with 187 additions and 222 deletions
ChangeLog
  View file @c1d193c
... ... @@ -6,6 +6,9 @@
6 6 Add equality and non-equality operators to GSkyRegionCircle class
7 7 Add GModelSpatialRadialRing class (#3253)
8 8 Fix ambiguity problem in GSkyDir::posang() and posang_deg() (#3199)
  9 + Fix GModelSpatialDiffuseCube pickeling issue (#3179)
  10 + Rename GModelSpatialDiffuseCube::set_mc_cone() to mc_cone() (#3179)
  11 + Rename GModelSpatialDiffuseMap::set_mc_cone() to mc_cone() (#3179)
9 12 Avoid rounding errors in GSkyMap operators (#3126)
10 13 Add GObservations::npred() methods (#3062)
11 14 Exclude upper edge boundary in GSkyMap::contains() method (#3012)
... ...
... ... @@ -48,6 +48,7 @@ The following methods have been added:
48 48 - GSkyRegionCircle::operator!=() - Non-equality operator
49 49 - GModelSpatial::flux() - Extract flux from sky region
50 50 - GModelSpatialDiffuseCube::mc_cone() - Return MC simulation cone
  51 +- GModelSpatialDiffuseMap::mc_cone() - Return MC simulation cone
51 52 - GCTAObservation::remove_response_cache(string&) - Remove response cache for model
52 53 - GCTAResponse::remove_response_cache(string&) - Remove response cache for model
53 54 - GCTAResponseCache::remove(string&) - Remove response cache for model
... ... @@ -62,6 +63,7 @@ The following methods have been removed:
62 63  
63 64 The following methods have been renamed:
64 65 - GModelSpatialDiffuseCube::set_mc_cone() to GModelSpatialDiffuseCube::mc_cone()
  66 +- GModelSpatialDiffuseMap::set_mc_cone() to GModelSpatialDiffuseMap::mc_cone()
65 67 - GResponse::irf(GEvent&, GSource&, GObservation&) to GResponse::irf_spatial()
66 68 - GLATResponse::irf(GLATEventAtom&, GSource&, GObservation&) to GLATResponse::irf_spatial_atom()
67 69 - GLATResponse::irf(GLATEventBin&, GSource&, GObservation&) to GLATResponse::irf_spatial_bin()
... ... @@ -110,6 +112,10 @@ None
110 112 The GModelSpatialRadialRing class was added that implements a radial ring model
111 113 with uniform source intensity (#3253).
112 114  
  115 +The set_mc_cone() methods of the GModelSpatialDiffuseCube and GModelSpatialDiffuseMap
  116 +class were renamed to mc_cone() and mc_cone() getter methods were implemented
  117 +for both classes. This solves a problem with pickling both classes (#3179).
  118 +
113 119 The method GModelSpatial::flux() was added that extracts the flux by integrating
114 120 over a specific sky region. For the moment only circular sky regions are
115 121 implemented (#2961).
... ...
doc/source/admin/release_history/1.7.rst
  View file @c1d193c
... ... @@ -19,7 +19,9 @@ Bug fixes
19 19 ---------
20 20  
21 21 * [`3199 <https://cta-redmine.irap.omp.eu/issues/3199>`_] -
22   - Fix ambiguity problem in ``GSkyDir::posang()`` and ``GSkyDir::posang_deg()``
  22 + Fix ambiguity problem in ``GSkyDir::posang()`` and ``GSkyDir::posang_deg()``.
  23 +* [`3179 <https://cta-redmine.irap.omp.eu/issues/3179>`_] -
  24 + Fix ``GModelSpatialDiffuseCube`` pickeling issue.
23 25 * [`3126 <https://cta-redmine.irap.omp.eu/issues/3126>`_] -
24 26 Avoid rounding errors in ``GSkyMap`` operators.
25 27 * [`3012 <https://cta-redmine.irap.omp.eu/issues/3012>`_] -
... ... @@ -82,5 +84,9 @@ New features
82 84 Clean up
83 85 --------
84 86  
  87 +* [`3179 <https://cta-redmine.irap.omp.eu/issues/3179>`_] -
  88 + Rename ``GModelSpatialDiffuseCube::set_mc_cone()`` to ``GModelSpatialDiffuseCube::mc_cone()`` and add ``mc_cone()``getter methods.
  89 +* [`3179 <https://cta-redmine.irap.omp.eu/issues/3179>`_] -
  90 + Rename ``GModelSpatialDiffuseMap::set_mc_cone()`` to ``GModelSpatialDiffuseMap::mc_cone()`` and add ``mc_cone()``getter methods.
85 91 * [`2893 <https://cta-redmine.irap.omp.eu/issues/2893>`_] -
86 92 Remove ``GCTAResponseCache::ndirs()``, ``GCTAResponseCache::nerecos()`` and ``GCTAResponseCache::netrues()`` methods.
... ...
include/GModelSpatialDiffuseCube.hpp
  View file @c1d193c
1 1 /***************************************************************************
2 2 * GModelSpatialDiffuseCube.hpp - Spatial map cube model class *
3 3 * ----------------------------------------------------------------------- *
4   - * copyright (C) 2010-2018 by Juergen Knoedlseder *
  4 + * copyright (C) 2010-2020 by Juergen Knoedlseder *
5 5 * ----------------------------------------------------------------------- *
6 6 * *
7 7 * This program is free software: you can redistribute it and/or modify *
... ... @@ -78,7 +78,6 @@ public:
78 78 virtual void clear(void);
79 79 virtual GModelSpatialDiffuseCube* clone(void) const;
80 80 virtual std::string classname(void) const;
81   - virtual std::string type(void) const;
82 81 virtual double eval(const GPhoton& photon,
83 82 const bool& gradients = false) const;
84 83 virtual GSkyDir mc(const GEnergy& energy,
... ... @@ -88,7 +87,6 @@ public:
88 87 const double& radius) const;
89 88 virtual bool contains(const GSkyDir& dir,
90 89 const double& margin = 0.0) const;
91   - virtual GSkyRegion* region(void) const;
92 90 virtual void read(const GXmlElement& xml);
93 91 virtual void write(GXmlElement& xml) const;
94 92 virtual std::string print(const GChatter& chatter = NORMAL) const;
... ... @@ -105,8 +103,8 @@ public:
105 103 GEnergies energies(void);
106 104 void energies(const GEnergies& energies);
107 105 const GModelSpectralNodes& spectrum(void) const;
108   - void set_mc_cone(const GSkyDir& centre,
109   - const double& radius) const;
  106 + void mc_cone(const GSkyRegionCircle& cone) const;
  107 + const GSkyRegionCircle& mc_cone(void) const;
110 108 void load(const GFilename& filename);
111 109 void save(const GFilename& filename,
112 110 const bool& clobber = false) const;
... ... @@ -115,29 +113,26 @@ public:
115 113  
116 114 protected:
117 115 // Protected methods
118   - void init_members(void);
119   - void copy_members(const GModelSpatialDiffuseCube& model);
120   - void free_members(void);
121   - void fetch_cube(void) const;
122   - void load_cube(const GFilename& filename);
123   - void set_energy_boundaries(void);
124   - void update_mc_cache(void);
125   - double cube_intensity(const GPhoton& photon) const;
126   - void set_region(void) const;
  116 + void init_members(void);
  117 + void copy_members(const GModelSpatialDiffuseCube& model);
  118 + void free_members(void);
  119 + void fetch_cube(void) const;
  120 + void load_cube(const GFilename& filename);
  121 + void set_energy_boundaries(void);
  122 + void update_mc_cache(void);
  123 + double cube_intensity(const GPhoton& photon) const;
  124 + virtual void set_region(void) const;
127 125  
128 126 // Protected members
129   - std::string m_type; //!< Model type
130   - GModelPar m_value; //!< Value
131   - GFilename m_filename; //!< Name of map cube
132   - bool m_loaded; //!< Signals if map cube has been loaded
133   - GSkyMap m_cube; //!< Map cube
134   - GNodeArray m_logE; //!< Log10(energy) values of the maps
135   - GEbounds m_ebounds; //!< Energy bounds of the maps
136   - mutable GSkyRegionCircle m_region; //!< Bounding circle
  127 + GModelPar m_value; //!< Value
  128 + GFilename m_filename; //!< Name of map cube
  129 + bool m_loaded; //!< Signals if map cube has been loaded
  130 + GSkyMap m_cube; //!< Map cube
  131 + GNodeArray m_logE; //!< Log10(energy) values of the maps
  132 + GEbounds m_ebounds; //!< Energy bounds of the maps
137 133  
138 134 // Monte Carlo cache
139   - mutable GSkyDir m_mc_centre; //!< Centre of MC cone
140   - mutable double m_mc_radius; //!< Radius of MC cone (degrees)
  135 + mutable GSkyRegionCircle m_mc_cone; //!< MC cone
141 136 mutable double m_mc_one_minus_cosrad; //!< 1-cosine of radius
142 137 mutable std::vector<double> m_mc_max; //!< Maximum values for MC
143 138 mutable GModelSpectralNodes m_mc_spectrum; //!< Map cube spectrum
... ... @@ -157,20 +152,6 @@ std::string GModelSpatialDiffuseCube::classname(void) const
157 152  
158 153  
159 154 /***********************************************************************//**
160   - * @brief Return spatial model type
161   - *
162   - * @return Model type.
163   - *
164   - * Returns the type of the spatial map cube model.
165   - ***************************************************************************/
166   -inline
167   -std::string GModelSpatialDiffuseCube::type(void) const
168   -{
169   - return (m_type);
170   -}
171   -
172   -
173   -/***********************************************************************//**
174 155 * @brief Return number of maps in cube
175 156 *
176 157 * @return Number of maps in cube.
... ... @@ -303,17 +284,16 @@ double GModelSpatialDiffuseCube::mc_norm(const GSkyDir&amp; dir,
303 284  
304 285  
305 286 /***********************************************************************//**
306   - * @brief Return boundary sky region
  287 + * @brief Get Monte Carlo simulation cone
307 288 *
308   - * @return Boundary sky region.
  289 + * @return Monte Carlo simulation sky region circle.
309 290 *
310   - * Returns a sky region that fully encloses the point source.
  291 + * Returns the sky region circle used for Monte Carlo simulations.
311 292 ***************************************************************************/
312 293 inline
313   -GSkyRegion* GModelSpatialDiffuseCube::region(void) const
  294 +const GSkyRegionCircle& GModelSpatialDiffuseCube::mc_cone(void) const
314 295 {
315   - set_region();
316   - return (&m_region);
  296 + return (m_mc_cone);
317 297 }
318 298  
319 299 #endif /* GMODELSPATIALDIFFUSECUBE_HPP */
... ...
include/GModelSpatialDiffuseMap.hpp
  View file @c1d193c
1 1 /***************************************************************************
2 2 * GModelSpatialDiffuseMap.hpp - Spatial map model class *
3 3 * ----------------------------------------------------------------------- *
4   - * copyright (C) 2012-2016 by Juergen Knoedlseder *
  4 + * copyright (C) 2012-2020 by Juergen Knoedlseder *
5 5 * ----------------------------------------------------------------------- *
6 6 * *
7 7 * This program is free software: you can redistribute it and/or modify *
... ... @@ -70,7 +70,6 @@ public:
70 70 virtual void clear(void);
71 71 virtual GModelSpatialDiffuseMap* clone(void) const;
72 72 virtual std::string classname(void) const;
73   - virtual std::string type(void) const;
74 73 virtual double eval(const GPhoton& photon,
75 74 const bool& gradients = false) const;
76 75 virtual GSkyDir mc(const GEnergy& energy,
... ... @@ -80,47 +79,41 @@ public:
80 79 const double& radius) const;
81 80 virtual bool contains(const GSkyDir& dir,
82 81 const double& margin = 0.0) const;
83   - virtual GSkyRegion* region(void) const;
84 82 virtual void read(const GXmlElement& xml);
85 83 virtual void write(GXmlElement& xml) const;
86 84 virtual std::string print(const GChatter& chatter = NORMAL) const;
87 85  
88 86 // Other methods
89   - double value(void) const;
90   - void value(const double& value);
91   - const GFilename& filename(void) const;
92   - void load(const GFilename& filename);
93   - const GSkyMap& map(void) const;
94   - void map(const GSkyMap& map);
95   - bool normalize(void) const;
96   - void set_mc_cone(const GSkyDir& centre,
97   - const double& radius) const;
  87 + double value(void) const;
  88 + void value(const double& value);
  89 + const GFilename& filename(void) const;
  90 + void load(const GFilename& filename);
  91 + const GSkyMap& map(void) const;
  92 + void map(const GSkyMap& map);
  93 + bool normalize(void) const;
  94 + void mc_cone(const GSkyRegionCircle& cone) const;
  95 + const GSkyRegionCircle& mc_cone(void) const;
98 96  
99 97 protected:
100 98 // Protected methods
101   - void init_members(void);
102   - void copy_members(const GModelSpatialDiffuseMap& model);
103   - void free_members(void);
104   - void prepare_map(void);
105   - void set_region(void) const;
  99 + void init_members(void);
  100 + void copy_members(const GModelSpatialDiffuseMap& model);
  101 + void free_members(void);
  102 + void prepare_map(void);
  103 + virtual void set_region(void) const;
106 104  
107 105 // Protected members
108   - std::string m_type; //!< Model type
109   - GModelPar m_value; //!< Value
110   - GSkyMap m_map; //!< Skymap
111   - GFilename m_filename; //!< Name of skymap
112   - GSkyDir m_centre; //!< Centre of bounding circle
113   - double m_radius; //!< Radius of bounding circle
114   - bool m_normalize; //!< Normalize map (default: true)
115   - bool m_has_normalize; //!< XML has normalize attribute
116   - mutable GSkyRegionCircle m_region; //!< Bounding circle
  106 + GModelPar m_value; //!< Value
  107 + GSkyMap m_map; //!< Skymap
  108 + GFilename m_filename; //!< Name of skymap
  109 + bool m_normalize; //!< Normalize map (default: true)
  110 + bool m_has_normalize; //!< XML has normalize attribute
117 111  
118 112 // MC simulation cache
119   - mutable GSkyDir m_mc_centre; //!< Centre of MC cone
120   - mutable double m_mc_radius; //!< Radius of MC cone
121   - mutable double m_mc_one_minus_cosrad; //!< 1-cosine of radius
122   - mutable double m_mc_norm; //!< Map normalization
123   - mutable double m_mc_max; //!< Map maximum for MC
  113 + mutable GSkyRegionCircle m_mc_cone; //!< MC cone
  114 + mutable double m_mc_one_minus_cosrad; //!< 1-cosine of radius
  115 + mutable double m_mc_norm; //!< Map normalization
  116 + mutable double m_mc_max; //!< Map maximum for MC
124 117 };
125 118  
126 119  
... ... @@ -137,20 +130,6 @@ std::string GModelSpatialDiffuseMap::classname(void) const
137 130  
138 131  
139 132 /***********************************************************************//**
140   - * @brief Return spatial model type
141   - *
142   - * @return Model type.
143   - *
144   - * Returns the type of the spatial map model.
145   - ***************************************************************************/
146   -inline
147   -std::string GModelSpatialDiffuseMap::type(void) const
148   -{
149   - return (m_type);
150   -}
151   -
152   -
153   -/***********************************************************************//**
154 133 * @brief Get model value
155 134 *
156 135 * @return Model value.
... ... @@ -238,17 +217,16 @@ bool GModelSpatialDiffuseMap::normalize(void) const
238 217  
239 218  
240 219 /***********************************************************************//**
241   - * @brief Return boundary sky region
  220 + * @brief Get Monte Carlo simulation cone
242 221 *
243   - * @return Boundary sky region.
  222 + * @return Monte Carlo simulation sky region circle.
244 223 *
245   - * Returns a sky region that fully encloses the point source.
  224 + * Returns the sky region circle used for Monte Carlo simulations.
246 225 ***************************************************************************/
247 226 inline
248   -GSkyRegion* GModelSpatialDiffuseMap::region(void) const
  227 +const GSkyRegionCircle& GModelSpatialDiffuseMap::mc_cone(void) const
249 228 {
250   - set_region();
251   - return (&m_region);
  229 + return (m_mc_cone);
252 230 }
253 231  
254 232 #endif /* GMODELSPATIALDIFFUSEMAP_HPP */
... ...
pyext/GModelSpatialDiffuseCube.i
  View file @c1d193c
1 1 /***************************************************************************
2 2 * GModelSpatialDiffuseCube.i - Spatial map cube model class *
3 3 * ----------------------------------------------------------------------- *
4   - * copyright (C) 2011-2018 by Juergen Knoedlseder *
  4 + * copyright (C) 2011-2020 by Juergen Knoedlseder *
5 5 * ----------------------------------------------------------------------- *
6 6 * *
7 7 * This program is free software: you can redistribute it and/or modify *
... ... @@ -52,7 +52,6 @@ public:
52 52 virtual void clear(void);
53 53 virtual GModelSpatialDiffuseCube* clone(void) const;
54 54 virtual std::string classname(void) const;
55   - virtual std::string type(void) const;
56 55 virtual double eval(const GPhoton& photon,
57 56 const bool& gradients = false) const;
58 57 virtual GSkyDir mc(const GEnergy& energy,
... ... @@ -62,7 +61,6 @@ public:
62 61 const double& radius) const;
63 62 virtual bool contains(const GSkyDir& dir,
64 63 const double& margin = 0.0) const;
65   - virtual GSkyRegion* region(void) const;
66 64 virtual void read(const GXmlElement& xml);
67 65 virtual void write(GXmlElement& xml) const;
68 66  
... ... @@ -78,8 +76,8 @@ public:
78 76 GEnergies energies(void);
79 77 void energies(const GEnergies& energies);
80 78 const GModelSpectralNodes& spectrum(void) const;
81   - void set_mc_cone(const GSkyDir& centre,
82   - const double& radius) const;
  79 + void mc_cone(const GSkyRegionCircle& cone) const;
  80 + const GSkyRegionCircle& mc_cone(void) const;
83 81 void load(const GFilename& filename);
84 82 void save(const GFilename& filename,
85 83 const bool& clobber = false) const;
... ... @@ -99,12 +97,13 @@ public:
99 97 def __getstate__(self):
100 98 fits = gammalib.GFits()
101 99 self.write(fits)
102   - state = (fits, self.filename(), self[0])
  100 + state = (fits, self.filename(), self[0], self.mc_cone())
103 101 return state
104 102 def __setstate__(self, state):
105 103 self.__init__()
106 104 self.read(state[0])
107 105 self.filename(state[1])
108 106 self[0] = state[2]
  107 + self.mc_cone(state[3])
109 108 }
110 109 };
... ...
pyext/GModelSpatialDiffuseMap.i
  View file @c1d193c
1 1 /***************************************************************************
2 2 * GModelSpatialDiffuseMap.i - Spatial map model class *
3 3 * ----------------------------------------------------------------------- *
4   - * copyright (C) 2012-2018 by Juergen Knoedlseder *
  4 + * copyright (C) 2012-2020 by Juergen Knoedlseder *
5 5 * ----------------------------------------------------------------------- *
6 6 * *
7 7 * This program is free software: you can redistribute it and/or modify *
... ... @@ -52,7 +52,6 @@ public:
52 52 virtual void clear(void);
53 53 virtual GModelSpatialDiffuseMap* clone(void) const;
54 54 virtual std::string classname(void) const;
55   - virtual std::string type(void) const;
56 55 virtual double eval(const GPhoton& photon,
57 56 const bool& gradients = false) const;
58 57 virtual GSkyDir mc(const GEnergy& energy,
... ... @@ -62,20 +61,19 @@ public:
62 61 const double& radius) const;
63 62 virtual bool contains(const GSkyDir& dir,
64 63 const double& margin = 0.0) const;
65   - virtual GSkyRegion* region(void) const;
66 64 virtual void read(const GXmlElement& xml);
67 65 virtual void write(GXmlElement& xml) const;
68 66  
69 67 // Other methods
70   - double value(void) const;
71   - void value(const double& value);
72   - const GFilename& filename(void) const;
73   - void load(const GFilename& filename);
74   - const GSkyMap& map(void) const;
75   - void map(const GSkyMap& map);
76   - bool normalize(void) const;
77   - void set_mc_cone(const GSkyDir& centre,
78   - const double& radius) const;
  68 + double value(void) const;
  69 + void value(const double& value);
  70 + const GFilename& filename(void) const;
  71 + void load(const GFilename& filename);
  72 + const GSkyMap& map(void) const;
  73 + void map(const GSkyMap& map);
  74 + bool normalize(void) const;
  75 + void mc_cone(const GSkyRegionCircle& cone) const;
  76 + const GSkyRegionCircle& mc_cone(void) const;
79 77 };
80 78  
81 79  
... ... @@ -90,9 +88,10 @@ public:
90 88 def __getstate__(self):
91 89 xml = gammalib.GXmlElement()
92 90 self.write(xml)
93   - state = (xml,)
  91 + state = (xml, self.mc_cone())
94 92 return state
95 93 def __setstate__(self, state):
96 94 self.__init__(state[0])
  95 + self.mc_cone(state[1])
97 96 }
98 97 };
... ...
src/model/GModelSky.cpp
  View file @c1d193c
1 1 /***************************************************************************
2 2 * GModelSky.cpp - Sky model class *
3 3 * ----------------------------------------------------------------------- *
4   - * copyright (C) 2011-2019 by Juergen Knoedlseder *
  4 + * copyright (C) 2011-2020 by Juergen Knoedlseder *
5 5 * ----------------------------------------------------------------------- *
6 6 * *
7 7 * This program is free software: you can redistribute it and/or modify *
... ... @@ -878,7 +878,7 @@ GPhotons GModelSky::mc(const double&amp; area,
878 878 if (cube != NULL) {
879 879  
880 880 // Set MC cone
881   - cube->set_mc_cone(dir, radius);
  881 + cube->mc_cone(GSkyRegionCircle(dir, radius));
882 882  
883 883 // Allocate node function to replace the spectral component
884 884 GModelSpectralNodes* nodes =
... ...
src/model/GModelSpatialDiffuseCube.cpp
  View file @c1d193c
1 1 /***************************************************************************
2 2 * GModelSpatialDiffuseCube.cpp - Spatial map cube model class *
3 3 * ----------------------------------------------------------------------- *
4   - * copyright (C) 2010-2018 by Juergen Knoedlseder *
  4 + * copyright (C) 2010-2020 by Juergen Knoedlseder *
5 5 * ----------------------------------------------------------------------- *
6 6 * *
7 7 * This program is free software: you can redistribute it and/or modify *
... ... @@ -458,7 +458,7 @@ GSkyDir GModelSpatialDiffuseCube::mc(const GEnergy&amp; energy,
458 458 double theta = std::acos(1.0 - ran.uniform() * m_mc_one_minus_cosrad) *
459 459 gammalib::rad2deg;
460 460 double phi = 360.0 * ran.uniform();
461   - dir = m_mc_centre;
  461 + dir = m_mc_cone.centre();
462 462 dir.rotate_deg(phi, theta);
463 463  
464 464 // Get map value at simulated sky direction
... ... @@ -761,26 +761,23 @@ GEnergies GModelSpatialDiffuseCube::energies(void)
761 761 /***********************************************************************//**
762 762 * @brief Set Monte Carlo simulation cone
763 763 *
764   - * @param[in] centre Simulation cone centre.
765   - * @param[in] radius Simulation cone radius (degrees).
  764 + * @param[in] cone Monte Carlo simulation cone.
766 765 *
767   - * Sets the simulation cone centre and radius that defines the directions
768   - * that will be simulated using the mc() method and pre-computes the maximum
769   - * intensity and the spatially integrated flux of each map within the
770   - * simulation cone region.
  766 + * Sets the simulation cone that defines the directions that will be
  767 + * simulated using the mc() method and pre-computes the maximum intensity
  768 + * and the spatially integrated flux of each map within the simulation cone
  769 + * region.
771 770 ***************************************************************************/
772   -void GModelSpatialDiffuseCube::set_mc_cone(const GSkyDir& centre,
773   - const double& radius) const
  771 +void GModelSpatialDiffuseCube::mc_cone(const GSkyRegionCircle& cone) const
774 772 {
775 773 // Continue only if the simulation cone has changed
776   - if ((centre != m_mc_centre) || (radius != m_mc_radius)) {
  774 + if (cone != m_mc_cone) {
777 775  
778 776 // Save simulation cone definition
779   - m_mc_centre = centre;
780   - m_mc_radius = radius;
  777 + m_mc_cone = cone;
781 778  
782 779 // Pre-compute 1 - cosine of radius
783   - m_mc_one_minus_cosrad = 1.0 - std::cos(m_mc_radius*gammalib::deg2rad);
  780 + m_mc_one_minus_cosrad = 1.0 - std::cos(m_mc_cone.radius()*gammalib::deg2rad);
784 781  
785 782 // Initialise Monte Carlo cache
786 783 m_mc_max.clear();
... ... @@ -804,8 +801,8 @@ void GModelSpatialDiffuseCube::set_mc_cone(const GSkyDir&amp; centre,
804 801 for (int k = 0; k < npix; ++k) {
805 802  
806 803 // Continue only if pixel is within MC cone
807   - double distance = centre.dist_deg(m_cube.pix2dir(k));
808   - if (distance <= radius) {
  804 + double distance = m_mc_cone.centre().dist_deg(m_cube.pix2dir(k));
  805 + if (distance <= m_mc_cone.radius()) {
809 806  
810 807 // Loop over all maps
811 808 for (int i = 0; i < nmaps; ++i) {
... ... @@ -857,7 +854,7 @@ void GModelSpatialDiffuseCube::set_mc_cone(const GSkyDir&amp; centre,
857 854  
858 855 // Log maximum intensity and total flux for debugging
859 856 #if defined(G_DEBUG_MC_CACHE)
860   - std::cout << "GModelSpatialDiffuseCube::set_mc_cone:" << std::endl;
  857 + std::cout << "GModelSpatialDiffuseCube::mc_cone:" << std::endl;
861 858 std::cout << " Maximum map intensity:" << std::endl;
862 859 for (int i = 0; i < m_mc_max.size(); ++i) {
863 860 GEnergy energy;
... ... @@ -1139,11 +1136,9 @@ void GModelSpatialDiffuseCube::init_members(void)
1139 1136 m_logE.clear();
1140 1137 m_ebounds.clear();
1141 1138 m_loaded = false;
1142   - m_region.clear();
1143 1139  
1144 1140 // Initialise MC cache
1145   - m_mc_centre.clear();
1146   - m_mc_radius = -1.0; //!< Signal that initialisation is needed
  1141 + m_mc_cone.clear();
1147 1142 m_mc_one_minus_cosrad = 1.0;
1148 1143 m_mc_max.clear();
1149 1144 m_mc_spectrum.clear();
... ... @@ -1161,18 +1156,16 @@ void GModelSpatialDiffuseCube::init_members(void)
1161 1156 void GModelSpatialDiffuseCube::copy_members(const GModelSpatialDiffuseCube& model)
1162 1157 {
1163 1158 // Copy members
1164   - m_type = model.m_type;
  1159 + m_type = model.m_type; // Needed to conserve model type
1165 1160 m_value = model.m_value;
1166 1161 m_filename = model.m_filename;
1167 1162 m_cube = model.m_cube;
1168 1163 m_logE = model.m_logE;
1169 1164 m_ebounds = model.m_ebounds;
1170 1165 m_loaded = model.m_loaded;
1171   - m_region = model.m_region;
1172 1166  
1173 1167 // Copy MC cache
1174   - m_mc_centre = model.m_mc_centre;
1175   - m_mc_radius = model.m_mc_radius;
  1168 + m_mc_cone = model.m_mc_cone;
1176 1169 m_mc_one_minus_cosrad = model.m_mc_one_minus_cosrad;
1177 1170 m_mc_max = model.m_mc_max;
1178 1171 m_mc_spectrum = model.m_mc_spectrum;
... ... @@ -1326,18 +1319,16 @@ void GModelSpatialDiffuseCube::set_energy_boundaries(void)
1326 1319 /***********************************************************************//**
1327 1320 * @brief Update Monte Carlo cache
1328 1321 *
1329   - * Initialise the cache for Monte Carlo sampling of the map cube. The Monte
1330   - * Carlo cache consists of a linear array that maps a value between 0 and 1
1331   - * into the skymap pixel for all maps in the cube.
  1322 + * Initialise the cache for Monte Carlo sampling of the map cube. See the
  1323 + * mc_cone() for more information.
1332 1324 ***************************************************************************/
1333 1325 void GModelSpatialDiffuseCube::update_mc_cache(void)
1334 1326 {
1335   - // Set centre and radius to all sky
1336   - GSkyDir centre;
1337   - double radius = 360.0;
  1327 + // Initialise sky region to all sky (dummy)
  1328 + GSkyRegionCircle cone(0.0, 0.0, 360.0);
1338 1329  
1339   - // Compute cache
1340   - set_mc_cone(centre, radius);
  1330 + // Set cone
  1331 + mc_cone(cone);
1341 1332  
1342 1333 // Return
1343 1334 return;
... ... @@ -1395,11 +1386,11 @@ double GModelSpatialDiffuseCube::cube_intensity(const GPhoton&amp; photon) const
1395 1386 ***************************************************************************/
1396 1387 void GModelSpatialDiffuseCube::set_region(void) const
1397 1388 {
1398   - // Set sky region centre to (0,0)
1399   - m_region.centre(0.0, 0.0);
  1389 + // Set sky region circle (all sky)
  1390 + GSkyRegionCircle region(0.0, 0.0, 180.0);
1400 1391  
1401   - // Set sky region radius to 180 degrees (all points included)
1402   - m_region.radius(180.0);
  1392 + // Set region (circumvent const correctness)
  1393 + const_cast<GModelSpatialDiffuseCube*>(this)->m_region = region;
1403 1394  
1404 1395 // Return
1405 1396 return;
... ...
src/model/GModelSpatialDiffuseMap.cpp
  View file @c1d193c
1 1 /***************************************************************************
2 2 * GModelSpatialDiffuseMap.cpp - Spatial map model class *
3 3 * ----------------------------------------------------------------------- *
4   - * copyright (C) 2012-2016 by Juergen Knoedlseder *
  4 + * copyright (C) 2012-2020 by Juergen Knoedlseder *
5 5 * ----------------------------------------------------------------------- *
6 6 * *
7 7 * This program is free software: you can redistribute it and/or modify *
... ... @@ -403,7 +403,7 @@ GSkyDir GModelSpatialDiffuseMap::mc(const GEnergy&amp; energy,
403 403 double theta = std::acos(1.0 - ran.uniform() * m_mc_one_minus_cosrad) *
404 404 gammalib::rad2deg;
405 405 double phi = 360.0 * ran.uniform();
406   - dir = m_mc_centre;
  406 + dir = m_mc_cone.centre();
407 407 dir.rotate_deg(phi, theta);
408 408  
409 409 // Get map value at simulated sky direction. If the map value is non-
... ... @@ -448,7 +448,7 @@ double GModelSpatialDiffuseMap::mc_norm(const GSkyDir&amp; dir,
448 448 const double& radius) const
449 449 {
450 450 // Set the MC cone
451   - set_mc_cone(dir, radius);
  451 + mc_cone(GSkyRegionCircle(dir, radius));
452 452  
453 453 // Retrieve normalization
454 454 double norm = m_mc_norm * value();
... ... @@ -475,14 +475,14 @@ bool GModelSpatialDiffuseMap::contains(const GSkyDir&amp; dir,
475 475 bool contains = false;
476 476  
477 477 // Continue only if radius is positive
478   - if (m_radius > 0.0) {
  478 + if (m_region.radius() > 0.0) {
479 479  
480 480 // Compute distance to centre
481   - double distance = m_centre.dist_deg(dir);
  481 + double distance = m_region.centre().dist_deg(dir);
482 482  
483 483 // If distance is smaller than radius plus margin we consider
484 484 // the position to be contained within the bounding circle
485   - if (distance < m_radius + margin) {
  485 + if (distance < m_region.radius() + margin) {
486 486 contains = true;
487 487 }
488 488  
... ... @@ -656,26 +656,23 @@ void GModelSpatialDiffuseMap::write(GXmlElement&amp; xml) const
656 656 /***********************************************************************//**
657 657 * @brief Set Monte Carlo simulation cone
658 658 *
659   - * @param[in] centre Simulation cone centre.
660   - * @param[in] radius Simulation cone radius (degrees).
  659 + * @param[in] cone Monte Carlo simulation cone.
661 660 *
662   - * Sets the simulation cone centre and radius that defines the directions
663   - * that will be simulated using the mc() method and pre-computes the maximum
664   - * intensity and the spatially integrated flux of the map within the
665   - * simulation cone region.
  661 + * Sets the simulation cone that defines the directions that will be
  662 + * simulated using the mc() method and pre-computes the maximum intensity and
  663 + * the spatially integrated flux of the map within the simulation cone
  664 + * region.
666 665 ***************************************************************************/
667   -void GModelSpatialDiffuseMap::set_mc_cone(const GSkyDir& centre,
668   - const double& radius) const
  666 +void GModelSpatialDiffuseMap::mc_cone(const GSkyRegionCircle& cone) const
669 667 {
670 668 // Continue only if the simulation cone has changed
671   - if ((centre != m_mc_centre) || (radius != m_mc_radius)) {
  669 + if (cone != m_mc_cone) {
672 670  
673 671 // Save simulation cone definition
674   - m_mc_centre = centre;
675   - m_mc_radius = radius;
  672 + m_mc_cone = cone;
676 673  
677 674 // Pre-compute 1 - cosine of radius
678   - m_mc_one_minus_cosrad = 1.0 - std::cos(m_mc_radius*gammalib::deg2rad);
  675 + m_mc_one_minus_cosrad = 1.0 - std::cos(m_mc_cone.radius()*gammalib::deg2rad);
679 676  
680 677 // Initialise map maximum and normalisation
681 678 m_mc_max = 0.0;
... ... @@ -695,18 +692,18 @@ void GModelSpatialDiffuseMap::set_mc_cone(const GSkyDir&amp; centre,
695 692 // Get map flux, intensity and distance from MC cone centre
696 693 double flux = m_map.flux(i);
697 694 double intensity = m_map(i);
698   - double distance = centre.dist_deg(m_map.pix2dir(i));
  695 + double distance = m_mc_cone.centre().dist_deg(m_map.pix2dir(i));
699 696  
700 697 // Add flux if positive
701 698 if (flux > 0.0) {
702   - if (distance <= radius) {
  699 + if (distance <= m_mc_cone.radius()) {
703 700 sum += flux; // flux within simulation cone
704 701 }
705 702 sum_map += flux; // total flux
706 703 }
707 704  
708 705 // Update maximum intensity
709   - if (distance <= radius) {
  706 + if (distance <= m_mc_cone.radius()) {
710 707 if (intensity > m_mc_max) {
711 708 m_mc_max = intensity;
712 709 }
... ... @@ -730,7 +727,7 @@ void GModelSpatialDiffuseMap::set_mc_cone(const GSkyDir&amp; centre,
730 727  
731 728 // Log maximum intensity and total flux for debugging
732 729 #if defined(G_DEBUG_MC_CACHE)
733   - std::cout << "GModelSpatialDiffuseMap::set_mc_cone:" << std::endl;
  730 + std::cout << "GModelSpatialDiffuseMap::mc_cone:" << std::endl;
734 731 std::cout << " Maximum map intensity:";
735 732 std::cout << m_mc_max << " ph/cm2/s/sr" << std::endl;
736 733 std::cout << " Spatially integrated flux:" << std::endl;
... ... @@ -774,9 +771,9 @@ std::string GModelSpatialDiffuseMap::print(const GChatter&amp; chatter) const
774 771 result.append(" [normalized]");
775 772 }
776 773 result.append("\n"+gammalib::parformat("Map centre"));
777   - result.append(m_centre.print());
  774 + result.append(m_region.centre().print());
778 775 result.append("\n"+gammalib::parformat("Map radius"));
779   - result.append(gammalib::str(m_radius)+" deg");
  776 + result.append(gammalib::str(m_region.radius())+" deg");
780 777 result.append("\n"+gammalib::parformat("Number of parameters"));
781 778 result.append(gammalib::str(size()));
782 779 for (int i = 0; i < size(); ++i) {
... ... @@ -860,13 +857,9 @@ void GModelSpatialDiffuseMap::init_members(void)
860 857 m_filename.clear();
861 858 m_normalize = true;
862 859 m_has_normalize = false;
863   - m_centre.clear();
864   - m_radius = 0.0;
865   - m_region.clear();
866 860  
867 861 // Initialise MC cache
868   - m_mc_centre.clear();
869   - m_mc_radius = -1.0; //!< Signal initialisation
  862 + m_mc_cone.clear();
870 863 m_mc_one_minus_cosrad = 1.0;
871 864 m_mc_norm = 0.0;
872 865 m_mc_max = 0.0;
... ... @@ -884,19 +877,15 @@ void GModelSpatialDiffuseMap::init_members(void)
884 877 void GModelSpatialDiffuseMap::copy_members(const GModelSpatialDiffuseMap& model)
885 878 {
886 879 // Copy members
887   - m_type = model.m_type;
  880 + m_type = model.m_type; // Needed to conserve model type
888 881 m_value = model.m_value;
889 882 m_map = model.m_map;
890 883 m_filename = model.m_filename;
891 884 m_normalize = model.m_normalize;
892 885 m_has_normalize = model.m_has_normalize;
893   - m_centre = model.m_centre;
894   - m_radius = model.m_radius;
895   - m_region = model.m_region;
896 886  
897 887 // Copy MC cache
898   - m_mc_centre = model.m_mc_centre;
899   - m_mc_radius = model.m_mc_radius;
  888 + m_mc_cone = model.m_mc_cone;
900 889 m_mc_one_minus_cosrad = model.m_mc_one_minus_cosrad;
901 890 m_mc_norm = model.m_mc_norm;
902 891 m_mc_max = model.m_mc_max;
... ... @@ -935,8 +924,8 @@ void GModelSpatialDiffuseMap::free_members(void)
935 924 void GModelSpatialDiffuseMap::prepare_map(void)
936 925 {
937 926 // Initialise centre and radius
938   - m_centre.clear();
939   - m_radius = 0.0;
  927 + m_region.clear();
  928 + //m_radius = 0.0;
940 929  
941 930 // Determine number of skymap pixels
942 931 int npix = m_map.npix();
... ... @@ -971,28 +960,34 @@ void GModelSpatialDiffuseMap::prepare_map(void)
971 960  
972 961 // If we have a HealPix map then set radius to 180 deg
973 962 if (m_map.projection()->code() == "HPX") {
974   - m_radius = 180.0;
  963 + m_region = GSkyRegionCircle(0.0, 0.0, 180.0);
975 964 }
976 965  
977 966 // ... otherwise compute map centre and radius
978 967 else {
979 968  
  969 + // Initialise maximum radius
  970 + double max_radius = 0.0;
  971 +
980 972 // Get map centre
981   - GSkyPixel centre(m_map.nx()/2.0, m_map.ny()/2.0);
982   - m_centre = m_map.pix2dir(centre);
  973 + GSkyPixel pixel(m_map.nx()/2.0, m_map.ny()/2.0);
  974 + GSkyDir centre = m_map.pix2dir(pixel);
983 975  
984 976 // Determine map radius
985 977 for (int i = 0; i < npix; ++i) {
986   - double radius = m_map.inx2dir(i).dist_deg(m_centre);
987   - if (radius > m_radius) {
988   - m_radius = radius;
  978 + double radius = m_map.inx2dir(i).dist_deg(centre);
  979 + if (radius > max_radius) {
  980 + max_radius = radius;
989 981 }
990 982 }
991 983  
  984 + // Set sky region
  985 + m_region = GSkyRegionCircle(centre, max_radius);
  986 +
992 987 } // endelse: computed map centre and radius
993 988  
994 989 // Set simulation cone
995   - set_mc_cone(m_centre, m_radius);
  990 + mc_cone(m_region);
996 991  
997 992 } // endif: there were skymap pixels
998 993  
... ... @@ -1006,11 +1001,13 @@ void GModelSpatialDiffuseMap::prepare_map(void)
1006 1001 ***************************************************************************/
1007 1002 void GModelSpatialDiffuseMap::set_region(void) const
1008 1003 {
1009   - // Set sky region centre to bounding circle centre
1010   - m_region.centre(m_centre);
  1004 + /*
  1005 + // Set sky region circle (all sky)
  1006 + GSkyRegionCircle region(m_centre, m_radius);
1011 1007  
1012   - // Set sky region radius to bounding circle radius
1013   - m_region.radius(m_radius);
  1008 + // Set region (circumvent const correctness)
  1009 + const_cast<GModelSpatialDiffuseMap*>(this)->m_region = region;
  1010 + */
1014 1011  
1015 1012 // Return
1016 1013 return;
... ...
test/test_GModel.py
  View file @c1d193c
... ... @@ -153,31 +153,37 @@ class Test(gammalib.GPythonTestSuite):
153 153 test_support.pickeling(self, gammalib.GModelTemporalPhaseCurve())
154 154  
155 155 # Setup for tests
156   - time = gammalib.GTime(1.0,'sec')
157   - pivot = gammalib.GEnergy(1.0,'TeV')
158   - ptsrc = gammalib.GModelSpatialPointSource(1.0,2.0)
159   - plaw = gammalib.GModelSpectralPlaw(1.0,-2.0,pivot)
160   - cnst = gammalib.GModelTemporalConst(3.0)
161   - sky = gammalib.GModelSky(ptsrc, plaw, cnst)
162   - engs = gammalib.GEnergies(10,pivot,gammalib.GEnergy(10.0,'TeV'))
163   - spacom = gammalib.GModelSpatialComposite()
  156 + fct = gammalib.GFilename(os.environ['TEST_DATA']+'/filefunction.txt')
  157 + lcrv = gammalib.GFilename(os.environ['TEST_DATA']+'/model_temporal_lightcurve.fits')
  158 + pcrv = gammalib.GFilename(os.environ['TEST_DATA']+'/model_temporal_phasecurve.fits')
  159 + cube = gammalib.GFilename(os.environ['TEST_DATA']+'/test_cube.fits')
  160 + table = gammalib.GFilename(os.environ['TEST_DATA']+'/model_point_table.fits')
  161 + dir = gammalib.GSkyDir()
  162 + dir.radec_deg(83.6331, 22.0145)
  163 + time = gammalib.GTime(1.0,'sec')
  164 + pivot = gammalib.GEnergy(1.0,'TeV')
  165 + ptsrc = gammalib.GModelSpatialPointSource(1.0,2.0)
  166 + plaw = gammalib.GModelSpectralPlaw(1.0,-2.0,pivot)
  167 + cnst = gammalib.GModelTemporalConst(3.0)
  168 + sky = gammalib.GModelSky(ptsrc, plaw, cnst)
  169 + engs = gammalib.GEnergies(10,pivot,gammalib.GEnergy(10.0,'TeV'))
  170 + circle = gammalib.GSkyRegionCircle(83.6331, 22.0145, 3.0)
  171 + spacom = gammalib.GModelSpatialComposite()
164 172 spacom.append(ptsrc, 'Src1')
165 173 spacom.append(ptsrc, 'Src2')
166   - speexp = gammalib.GModelSpectralExponential()
  174 + speexp = gammalib.GModelSpectralExponential()
167 175 speexp.exponent(plaw)
168   - specom = gammalib.GModelSpectralComposite()
  176 + specom = gammalib.GModelSpectralComposite()
169 177 specom.append(plaw, 'Src1')
170 178 specom.append(plaw, 'Src2')
171   - spemul = gammalib.GModelSpectralMultiplicative()
  179 + spemul = gammalib.GModelSpectralMultiplicative()
172 180 spemul.append(plaw, 'Src1')
173 181 spemul.append(plaw, 'Src2')
174   - nodes = gammalib.GModelSpectralNodes(plaw, engs)
175   - fct = gammalib.GFilename(os.environ['TEST_DATA']+'/filefunction.txt')
176   - lcrv = gammalib.GFilename(os.environ['TEST_DATA']+'/model_temporal_lightcurve.fits')
177   - pcrv = gammalib.GFilename(os.environ['TEST_DATA']+'/model_temporal_phasecurve.fits')
178   - cube = gammalib.GFilename(os.environ['TEST_DATA']+'/test_cube.fits')
179   - table = gammalib.GFilename(os.environ['TEST_DATA']+'/model_point_table.fits')
180   - dir = gammalib.GSkyDir()
  182 + spenod = gammalib.GModelSpectralNodes(plaw, engs)
  183 + spamap = gammalib.GModelSpatialDiffuseMap(cube, 2.0)
  184 + spamap.mc_cone(circle)
  185 + spacube = gammalib.GModelSpatialDiffuseCube(cube, 3.0)
  186 + spacube.mc_cone(circle)
181 187  
182 188 # Perform pickeling tests of filled classes
183 189 test_support.pickeling(self, gammalib.GModelPar('Par',2.0,5.0))
... ... @@ -185,8 +191,8 @@ class Test(gammalib.GPythonTestSuite):
185 191 test_support.pickeling(self, gammalib.GModelSky(sky))
186 192 test_support.pickeling(self, gammalib.GModelSpatialComposite(spacom))
187 193 test_support.pickeling(self, gammalib.GModelSpatialDiffuseConst(2.0))
188   - test_support.pickeling(self, gammalib.GModelSpatialDiffuseCube(cube, 2.0))
189   - test_support.pickeling(self, gammalib.GModelSpatialDiffuseMap(cube,2.0))
  194 + test_support.pickeling(self, spacube)
  195 + test_support.pickeling(self, spamap)
190 196 test_support.pickeling(self, gammalib.GModelSpatialEllipticalDisk(dir,0.2,0.3,45.0))
191 197 test_support.pickeling(self, gammalib.GModelSpatialEllipticalGauss(dir,0.2,0.3,45.0))
192 198 test_support.pickeling(self, gammalib.GModelSpatialPointSource(ptsrc))
... ... @@ -208,7 +214,7 @@ class Test(gammalib.GPythonTestSuite):
208 214 test_support.pickeling(self, gammalib.GModelSpectralGauss(1.0,pivot,pivot))
209 215 test_support.pickeling(self, gammalib.GModelSpectralLogParabola(1.0,-2.0,pivot,0.9))
210 216 test_support.pickeling(self, gammalib.GModelSpectralMultiplicative(spemul))
211   - test_support.pickeling(self, gammalib.GModelSpectralNodes(nodes))
  217 + test_support.pickeling(self, gammalib.GModelSpectralNodes(spenod))
212 218 test_support.pickeling(self, gammalib.GModelSpectralPlaw(plaw))
213 219 test_support.pickeling(self, gammalib.GModelSpectralPlawEnergyFlux(1.0,-2.0,pivot,pivot))
214 220 test_support.pickeling(self, gammalib.GModelSpectralPlawPhotonFlux(1.0,-2.0,pivot,pivot))
... ...