Eirik Hatlen / gammalib

Commit 89470207a17861e7ab60a0e0edad3a92eb89a09c

Authored by Jürgen Knödlseder
1 parent 600da647

Integrate code

Inline Side-by-side
Showing 14 changed files with 308 additions and 293 deletions
... ... @@ -15,4 +15,5 @@ Michael Mayer (michael.mayer@physik.hu-berlin.de)
15 15 Ellis Owen (Ellis.Owen@mpi-hd.mpg.de)
16 16 Anneli Schulz (anneli.schulz@desy.de)
17 17 Hubert Siejkowski (h.siejkowski@cyfronet.pl)
18   -Alexander Ziegler (alexander.ziegler@fau.de)
19 18 \ No newline at end of file
  19 +Domenico Tiziani (domenico.tiziani@fau.de)
  20 +Alexander Ziegler (alexander.ziegler@fau.de)
... ...
include/GModelSpatialComposite.hpp
  View file @8947020
... ... @@ -42,7 +42,6 @@
42 42 *
43 43 * This class implements the spatial model for a composition of spatial
44 44 * models.
45   - *
46 45 ***************************************************************************/
47 46 class GModelSpatialComposite : public GModelSpatial {
48 47  
... ... @@ -58,27 +57,28 @@ public:
58 57 virtual GModelSpatialComposite& operator=(const GModelSpatialComposite& model);
59 58  
60 59 // Implemented pure virtual methods
61   - virtual void clear(void);
  60 + virtual void clear(void);
62 61 virtual GModelSpatialComposite* clone(void) const;
63   - virtual std::string classname(void) const;
64   - virtual std::string type(void) const;
65   - virtual GClassCode code(void) const;
66   - virtual double eval(const GPhoton& photon,
67   - const bool& gradients = false) const;
68   - virtual GSkyDir mc(const GEnergy& energy,
69   - const GTime& time,
70   - GRan& ran) const;
71   - virtual double mc_norm(const GSkyDir& dir,
72   - const double& radius) const;
73   - virtual bool contains(const GSkyDir& dir,
74   - const double& margin = 0.0) const;
75   - virtual void read(const GXmlElement& xml);
76   - virtual void write(GXmlElement& xml) const;
77   - virtual std::string print(const GChatter& chatter = NORMAL) const;
  62 + virtual std::string classname(void) const;
  63 + virtual std::string type(void) const;
  64 + virtual GClassCode code(void) const;
  65 + virtual double eval(const GPhoton& photon,
  66 + const bool& gradients = false) const;
  67 + virtual GSkyDir mc(const GEnergy& energy,
  68 + const GTime& time,
  69 + GRan& ran) const;
  70 + virtual double mc_norm(const GSkyDir& dir,
  71 + const double& radius) const;
  72 + virtual bool contains(const GSkyDir& dir,
  73 + const double& margin = 0.0) const;
  74 + virtual void read(const GXmlElement& xml);
  75 + virtual void write(GXmlElement& xml) const;
  76 + virtual std::string print(const GChatter& chatter = NORMAL) const;
78 77  
79 78 // Other methods
80   - int components(void) const;
81   - void append(const GModelSpatial& component, const std::string& name = "");
  79 + int components(void) const;
  80 + void append(const GModelSpatial& component,
  81 + const std::string& name = "");
82 82 const GModelSpatial* component(const int& index) const;
83 83 const GModelSpatial* component(const std::string& name) const;
84 84  
... ... @@ -90,9 +90,9 @@ protected:
90 90 void free_members(void);
91 91  
92 92 // Protected members
93   - std::string m_type; //!< Model type
94   - std::vector<GModelSpatial*> m_components; //!< Components
95   - std::vector<std::string> m_component_names; //!< Component names
  93 + std::string m_type; //!< Model type
  94 + std::vector<GModelSpatial*> m_components; //!< Components
  95 + std::vector<std::string> m_component_names; //!< Component names
96 96 };
97 97  
98 98  
... ... @@ -135,24 +135,34 @@ GClassCode GModelSpatialComposite::code(void) const
135 135 return GMODEL_SPATIAL_COMPOSITE;
136 136 }
137 137  
  138 +
138 139 /***********************************************************************//**
139   - * @brief Return normalization of point source for Monte Carlo simulations
  140 + * @brief Return normalization of composite source for Monte Carlo
  141 + * simulations
140 142 *
141 143 * @param[in] dir Centre of simulation cone.
142 144 * @param[in] radius Radius of simulation cone (degrees).
143 145 * @return Normalization.
144 146 *
145   - * Returns the normalization for a point source within a circular region.
146   - * The normalization is 1 if the point source falls within the circle
147   - * defined by @p dir and @p radius, 0 otherwise.
  147 + * Returns unity.
148 148 ***************************************************************************/
149 149 inline
150 150 double GModelSpatialComposite::mc_norm(const GSkyDir& dir,
151   - const double& radius) const
  151 + const double& radius) const
152 152 {
153   - //double norm = (dir.dist_deg(this->dir()) <= radius) ? 1.0 : 0.0;
154   - //return (norm);
155 153 return 1.0;
156 154 }
157 155  
  156 +
  157 +/***********************************************************************//**
  158 + * @brief Return number of model components
  159 + *
  160 + * @return Number of model components.
  161 + ***************************************************************************/
  162 +inline
  163 +int GModelSpatialComposite::components(void) const
  164 +{
  165 + return m_components.size();
  166 +}
  167 +
158 168 #endif /* GMODELSPATIALCOMPOSITE_HPP */
... ...
include/GTools.hpp
  View file @8947020
... ... @@ -105,7 +105,7 @@ namespace gammalib {
105 105 bool contains(const std::string& str,
106 106 const std::string& substring);
107 107 bool contains(const std::vector<std::string> strings,
108   - const std::string& string);
  108 + const std::string& string);
109 109 void warning(const std::string& origin,
110 110 const std::string& message);
111 111 std::string xml2str(const std::string& arg);
... ...
inst/cta/include/GCTAResponseCube.hpp
  View file @8947020
... ... @@ -137,8 +137,8 @@ private:
137 137 const GSource& source,
138 138 const GObservation& obs) const;
139 139 double irf_composite(const GEvent& event,
140   - const GSource& source,
141   - const GObservation& obs) const;
  140 + const GSource& source,
  141 + const GObservation& obs) const;
142 142  
143 143 // Private data members
144 144 GCTACubeExposure m_exposure; //!< Exposure cube
... ...
inst/cta/include/GCTAResponseIrf.hpp
  View file @8947020
... ... @@ -177,8 +177,8 @@ private:
177 177 const GSource& source,
178 178 const GObservation& obs) const;
179 179 double irf_composite(const GEvent& event,
180   - const GSource& source,
181   - const GObservation& obs) const;
  180 + const GSource& source,
  181 + const GObservation& obs) const;
182 182 double nroi_ptsrc(const GModelSky& model,
183 183 const GEnergy& srcEng,
184 184 const GTime& srcTime,
... ... @@ -204,11 +204,11 @@ private:
204 204 const GTime& obsTime,
205 205 const GObservation& obs) const;
206 206 double nroi_composite(const GModelSky& model,
207   - const GEnergy& srcEng,
208   - const GTime& srcTime,
209   - const GEnergy& obsEng,
210   - const GTime& obsTime,
211   - const GObservation& obs) const;
  207 + const GEnergy& srcEng,
  208 + const GTime& srcTime,
  209 + const GEnergy& obsEng,
  210 + const GTime& obsTime,
  211 + const GObservation& obs) const;
212 212  
213 213 // Private data members
214 214 GCaldb m_caldb; //!< Calibration database
... ...
inst/cta/src/GCTAResponseCube.cpp
  View file @8947020
... ... @@ -1510,17 +1510,18 @@ double GCTAResponseCube::irf_diffuse(const GEvent&amp; event,
1510 1510 * Returns the instrument response to a specified composite source.
1511 1511 ***************************************************************************/
1512 1512 double GCTAResponseCube::irf_composite(const GEvent& event,
1513   - const GSource& source,
1514   - const GObservation& obs) const
  1513 + const GSource& source,
  1514 + const GObservation& obs) const
1515 1515 {
1516 1516 // Initialise IRF
1517 1517 double irf = 0.0;
1518 1518  
1519 1519 // Get pointer to composite model
1520   - const GModelSpatialComposite* model = dynamic_cast<const GModelSpatialComposite*>(source.model());
  1520 + const GModelSpatialComposite* model =
  1521 + dynamic_cast<const GModelSpatialComposite*>(source.model());
1521 1522  
1522 1523 // Loop over model components
1523   - for (int i=0;i<model->components();++i) {
  1524 + for (int i = 0; i < model->components(); ++i) {
1524 1525  
1525 1526 // Get pointer to spatial component
1526 1527 GModelSpatial* spat = const_cast<GModelSpatial*>(model->component(i));
... ... @@ -1530,11 +1531,13 @@ double GCTAResponseCube::irf_composite(const GEvent&amp; event,
1530 1531  
1531 1532 // Compute irf value
1532 1533 irf += this->irf(event, source, obs);
  1534 +
1533 1535 }
1534 1536  
1535 1537 // Divide by number of model components
1536   - irf /= (double)model->components();
1537   -
  1538 + if (model->components() > 0) {
  1539 + irf /= double(model->components());
  1540 + }
1538 1541  
1539 1542 return irf;
1540 1543 }
... ...
inst/cta/src/GCTAResponseIrf.cpp
  View file @8947020
... ... @@ -2863,14 +2863,15 @@ double GCTAResponseIrf::irf_diffuse(const GEvent&amp; event,
2863 2863 * Returns the instrument response to a specified composite source.
2864 2864 ***************************************************************************/
2865 2865 double GCTAResponseIrf::irf_composite(const GEvent& event,
2866   - const GSource& source,
2867   - const GObservation& obs) const
  2866 + const GSource& source,
  2867 + const GObservation& obs) const
2868 2868 {
2869 2869 // Initialise IRF
2870 2870 double irf = 0.0;
2871 2871  
2872 2872 // Get pointer to composite model
2873   - const GModelSpatialComposite* model = dynamic_cast<const GModelSpatialComposite*>(source.model());
  2873 + const GModelSpatialComposite* model =
  2874 + dynamic_cast<const GModelSpatialComposite*>(source.model());
2874 2875  
2875 2876 // Loop over model components
2876 2877 for (int i = 0; i < model->components(); ++i) {
... ... @@ -2883,10 +2884,13 @@ double GCTAResponseIrf::irf_composite(const GEvent&amp; event,
2883 2884  
2884 2885 // Compute irf value
2885 2886 irf += this->irf(event, src, obs);
  2887 +
2886 2888 }
2887 2889  
2888 2890 // Divide by number of model components
2889   - irf /= (double)model->components();
  2891 + if (model->components() > 0) {
  2892 + irf /= double(model->components());
  2893 + }
2890 2894  
2891 2895 // Return IRF value
2892 2896 return irf;
... ... @@ -3509,17 +3513,18 @@ double GCTAResponseIrf::nroi_diffuse(const GModelSky&amp; model,
3509 3513 * \f$S(p,E,t)\f$ and the response function \f$R(p',E',t'|p,E,t)\f$.
3510 3514 ***************************************************************************/
3511 3515 double GCTAResponseIrf::nroi_composite(const GModelSky& model,
3512   - const GEnergy& srcEng,
3513   - const GTime& srcTime,
3514   - const GEnergy& obsEng,
3515   - const GTime& obsTime,
3516   - const GObservation& obs) const
  3516 + const GEnergy& srcEng,
  3517 + const GTime& srcTime,
  3518 + const GEnergy& obsEng,
  3519 + const GTime& obsTime,
  3520 + const GObservation& obs) const
3517 3521 {
3518 3522 // Initialise nroi
3519 3523 double nroi = 0.0;
3520 3524  
3521 3525 // Get composite model
3522   - GModelSpatialComposite* comp = dynamic_cast<GModelSpatialComposite*>(model.spatial());
  3526 + GModelSpatialComposite* comp =
  3527 + dynamic_cast<GModelSpatialComposite*>(model.spatial());
3523 3528  
3524 3529 // Loop over model components
3525 3530 for (int i = 0; i < comp->components(); ++i) {
... ... @@ -3532,10 +3537,13 @@ double GCTAResponseIrf::nroi_composite(const GModelSky&amp; model,
3532 3537  
3533 3538 // Compute nroi
3534 3539 nroi += this->nroi(sky, srcEng, srcTime, obsEng, obsTime, obs);
  3540 +
3535 3541 }
3536 3542  
3537 3543 // Divide by number of model components
3538   - nroi /= (double)comp->components();
  3544 + if (comp->components() > 0) {
  3545 + nroi /= double(comp->components());
  3546 + }
3539 3547  
3540 3548 // Return nroi
3541 3549 return nroi;
... ...
pyext/GModelSpatialComposite.i
  View file @8947020
1 1 /***************************************************************************
2   - * GModelSpatialComposite.i - Spatial point source model class *
  2 + * GModelSpatialComposite.i - Spatial point source model class *
3 3 * ----------------------------------------------------------------------- *
4 4 * copyright (C) 2016 by Domenico Tiziani *
5 5 * ----------------------------------------------------------------------- *
... ... @@ -44,26 +44,26 @@ public:
44 44 virtual ~GModelSpatialComposite(void);
45 45  
46 46 // Implemented virtual methods
47   - virtual void clear(void);
  47 + virtual void clear(void);
48 48 virtual GModelSpatialComposite* clone(void) const;
49   - virtual std::string classname(void) const;
50   - virtual std::string type(void) const;
51   - virtual double eval(const GPhoton& photon,
52   - const bool& gradients = false) const;
53   - virtual GSkyDir mc(const GEnergy& energy,
54   - const GTime& time,
55   - GRan& ran) const;
56   - virtual double mc_norm(const GSkyDir& dir,
57   - const double& radius) const;
58   - virtual bool contains(const GSkyDir& dir,
59   - const double& margin = 0.0) const;
60   - virtual void read(const GXmlElement& xml);
61   - virtual void write(GXmlElement& xml) const;
  49 + virtual std::string classname(void) const;
  50 + virtual std::string type(void) const;
  51 + virtual double eval(const GPhoton& photon,
  52 + const bool& gradients = false) const;
  53 + virtual GSkyDir mc(const GEnergy& energy,
  54 + const GTime& time,
  55 + GRan& ran) const;
  56 + virtual double mc_norm(const GSkyDir& dir,
  57 + const double& radius) const;
  58 + virtual bool contains(const GSkyDir& dir,
  59 + const double& margin = 0.0) const;
  60 + virtual void read(const GXmlElement& xml);
  61 + virtual void write(GXmlElement& xml) const;
62 62  
63 63 // Other methods
64   - int components(void) const;
65   - void append(const GModelSpatial& component, const std::string& name = "");
66   - void append(const GModelSpatial& component);
  64 + int components(void) const;
  65 + void append(const GModelSpatial& component,
  66 + const std::string& name = "");
67 67 const GModelSpatial* component(const int& index) const;
68 68 const GModelSpatial* component(const std::string& name) const;
69 69 };
... ...
src/model/GModelSpatialComposite.cpp
  View file @8947020
1 1 /***************************************************************************
2   - * GModelSpatialComposite.cpp - Spatial point source model class *
  2 + * GModelSpatialComposite.cpp - Spatial composite model class *
3 3 * ----------------------------------------------------------------------- *
4 4 * copyright (C) 2016 by Domenico Tiziani *
5 5 * ----------------------------------------------------------------------- *
... ... @@ -39,14 +39,15 @@
39 39  
40 40 /* __ Globals ____________________________________________________________ */
41 41 const GModelSpatialComposite g_spatial_comp_seed;
42   -const GModelSpatialRegistry g_spatial_comp_registry(&g_spatial_comp_seed);
  42 +const GModelSpatialRegistry g_spatial_comp_registry(&g_spatial_comp_seed);
43 43  
44 44 /* __ Method name definitions ____________________________________________ */
45   -#define G_READ "GModelSpatialComposite::read(GXmlElement&)"
46   -#define G_WRITE "GModelSpatialComposite::write(GXmlElement&)"
47   -#define G_COMPONENT_INDEX "GModelSpatialComposite::component(const int&)"
48   -#define G_COMPONENT_NAME "GModelSpatialComposite::component(const std::string&)"
49   -#define G_APPEND "GModelSpatialComposite::append(const GModelSpatial&, const std::string&)"
  45 +#define G_READ "GModelSpatialComposite::read(GXmlElement&)"
  46 +#define G_WRITE "GModelSpatialComposite::write(GXmlElement&)"
  47 +#define G_COMPONENT_INDEX "GModelSpatialComposite::component(int&)"
  48 +#define G_COMPONENT_NAME "GModelSpatialComposite::component(std::string&)"
  49 +#define G_APPEND "GModelSpatialComposite::append(GModelSpatial&, "\
  50 + "std::string&)"
50 51  
51 52 /* __ Macros _____________________________________________________________ */
52 53  
... ... @@ -85,8 +86,8 @@ GModelSpatialComposite::GModelSpatialComposite(void) : GModelSpatial()
85 86 * Constructs empty spatial composite model by specifying a model @p type.
86 87 ***************************************************************************/
87 88 GModelSpatialComposite::GModelSpatialComposite(const bool& dummy,
88   - const std::string& type) :
89   - GModelSpatial()
  89 + const std::string& type) :
  90 + GModelSpatial()
90 91 {
91 92 // Initialise members
92 93 init_members();
... ... @@ -98,6 +99,7 @@ GModelSpatialComposite::GModelSpatialComposite(const bool&amp; dummy,
98 99 return;
99 100 }
100 101  
  102 +
101 103 /***********************************************************************//**
102 104 * @brief XML constructor
103 105 *
... ... @@ -108,7 +110,7 @@ GModelSpatialComposite::GModelSpatialComposite(const bool&amp; dummy,
108 110 * structure of the XML element.
109 111 ***************************************************************************/
110 112 GModelSpatialComposite::GModelSpatialComposite(const GXmlElement& xml) :
111   - GModelSpatial()
  113 + GModelSpatial()
112 114 {
113 115 // Initialise members
114 116 init_members();
... ... @@ -127,7 +129,7 @@ GModelSpatialComposite::GModelSpatialComposite(const GXmlElement&amp; xml) :
127 129 * @param[in] model Spatial composite model.
128 130 ***************************************************************************/
129 131 GModelSpatialComposite::GModelSpatialComposite(const GModelSpatialComposite& model) :
130   - GModelSpatial(model)
  132 + GModelSpatial(model)
131 133 {
132 134 // Initialise members
133 135 init_members();
... ... @@ -232,24 +234,24 @@ GModelSpatialComposite* GModelSpatialComposite::clone(void) const
232 234 * @param[in] gradients Compute gradients?
233 235 * @return Model value.
234 236 *
235   - * Evaluates the spatial part for all model components and returns
236   - * the normalised sum.
237   - *
  237 + * Evaluates the spatial composite model by summing all model components and
  238 + * dividing the sum by the number of components.
238 239 ***************************************************************************/
239 240 double GModelSpatialComposite::eval(const GPhoton& photon,
240   - const bool& gradients) const
  241 + const bool& gradients) const
241 242 {
242 243 // Initalise value
243   - double value = 0.;
  244 + double value = 0.0;
244 245  
245 246 // Sum over all components
246   - for (unsigned int i = 0; i < components(); i++) {
  247 + for (int i = 0; i < m_components.size(); ++i) {
247 248 value += m_components[i]->eval(photon, gradients);
248 249 }
249 250  
250   - // Normalise
251   - value /= (double)m_components.size();
252   -
  251 + // Normalise sum by the number of components
  252 + if (m_components.size() > 0) {
  253 + value /= double(m_components.size());
  254 + }
253 255  
254 256 // Return value
255 257 return value;
... ... @@ -268,24 +270,17 @@ double GModelSpatialComposite::eval(const GPhoton&amp; photon,
268 270 * component.
269 271 ***************************************************************************/
270 272 GSkyDir GModelSpatialComposite::mc(const GEnergy& energy,
271   - const GTime& time,
272   - GRan& ran) const
  273 + const GTime& time,
  274 + GRan& ran) const
273 275 {
274   - // Choose one component
275   -
276   - int idx = 1;
277   - double u = ran.uniform() * m_components.size();
278   -
279   - while (idx < m_components.size()) {
280   - if (u <= idx) {
281   - break;
282   - }
283   - idx++;
284   - }
285   - GModelSpatial *component = m_components[idx-1];
  276 + // Randomly choose one model component
  277 + int index = int(ran.uniform() * double(m_components.size()));
  278 + if (index >= m_components.size()) {
  279 + index = m_components.size();
  280 + }
286 281  
287   - // Simulate component
288   - GSkyDir sky_dir = component->mc(energy, time, ran);
  282 + // Draw random sky direction from the selected component
  283 + GSkyDir sky_dir = m_components[index]->mc(energy, time, ran);
289 284  
290 285 // Return sky direction
291 286 return (sky_dir);
... ... @@ -303,14 +298,21 @@ GSkyDir GModelSpatialComposite::mc(const GEnergy&amp; energy,
303 298 * model components.
304 299 ***************************************************************************/
305 300 bool GModelSpatialComposite::contains(const GSkyDir& dir,
306   - const double& margin) const
  301 + const double& margin) const
307 302 {
  303 + // Initialise containment flag
  304 + bool containment = false;
  305 +
308 306 // Loop over all components
309   - for (unsigned int idx_comp = 0; idx_comp < components(); idx_comp++) {
310   - if (m_components[idx_comp]->contains(dir, margin))
311   - return true;
  307 + for (int i = 0; i < m_components.size(); ++i) {
  308 + if (m_components[i]->contains(dir, margin)) {
  309 + containment = true;
  310 + break;
  311 + }
312 312 }
313   - return false;
  313 +
  314 + // Return containment
  315 + return containment;
314 316 }
315 317  
316 318  
... ... @@ -327,7 +329,7 @@ void GModelSpatialComposite::read(const GXmlElement&amp; xml)
327 329 int n_comp = xml.elements("spatialModel");
328 330  
329 331 // Loop over spatial elements
330   - for(int i = 0; i < n_comp; ++i) {
  332 + for (int i = 0; i < n_comp; ++i) {
331 333  
332 334 // Get spatial XML element
333 335 const GXmlElement* spec = xml.element("spatialModel", i);
... ... @@ -363,16 +365,7 @@ void GModelSpatialComposite::read(const GXmlElement&amp; xml)
363 365 * @exception GException::model_invalid_parnames
364 366 * Invalid model parameter names found in XML element.
365 367 *
366   - * Write the point source information into an XML element with the following
367   - * format
368   - *
369   - * <spatialModel type="PointSource">
370   - * <parameter free="0" max="360" min="-360" name="RA" scale="1" value="83.6331" />
371   - * <parameter free="0" max="90" min="-90" name="DEC" scale="1" value="22.0145" />
372   - * </spatialModel>
373   - *
374   - * @todo The case that an existing spatial XML element with "GLON" and "GLAT"
375   - * as coordinates is not supported.
  368 + * Write the spatial information into an XML element.
376 369 ***************************************************************************/
377 370 void GModelSpatialComposite::write(GXmlElement& xml) const
378 371 {
... ... @@ -388,65 +381,61 @@ void GModelSpatialComposite::write(GXmlElement&amp; xml) const
388 381 }
389 382  
390 383 // Write all model components
391   - for (unsigned int idx_comp = 0; idx_comp < m_components.size(); ++idx_comp) {
392   - GModelSpatial *component = m_components[idx_comp];
  384 + for (int i = 0; i < m_components.size(); ++i) {
  385 +
  386 + // Get spatial component
  387 + GModelSpatial* component = m_components[i];
393 388  
394   - // Check for NULL
395   - if (!component)
  389 + // Fall through if component is empty
  390 + if (component == NULL) {
396 391 continue;
  392 + }
397 393  
398   - // Find xml element with matching name
  394 + // Find XML element with matching name
399 395 GXmlElement *matching_model = NULL;
400   - for (unsigned int idx_elt = 0; idx_elt < xml.elements("spatialModel"); idx_elt++) {
401   - if (xml.element("spatialModel", idx_elt)->attribute("name") ==
402   - m_component_names[idx_comp]) {
403   - matching_model = xml.element("spatialModel", idx_elt);
  396 + for (int k = 0; k < xml.elements("spatialModel"); ++k) {
  397 + if (xml.element("spatialModel", k)->attribute("name") ==
  398 + m_component_names[i]) {
  399 + matching_model = xml.element("spatialModel", k);
404 400 break;
405 401 }
406   - } // endfor: loop over all xml elements
407   -
408   - // Create temporary copy of the spatial model
409   - GModelSpatial* spat = component->clone();
  402 + } // endfor: loop over all XML elements
410 403  
411   - // Strip prefix from parameter names
  404 + // Create copy of the spatial model
  405 + GModelSpatial* spatial = component->clone();
412 406  
413   - // Loop over all parameters of component
414   - for (unsigned int idx_par = 0; idx_par < spat->size(); ++idx_par) {
  407 + // Loop over all parameters of component and strip prefix
  408 + for (int k = 0; k < spatial->size(); ++k) {
415 409  
416   - // Get parameter
417   - GModelPar& par = (*spat)[idx_par];
418   -
419   - // Check if name contains colon
420   - if (gammalib::contains(par.name(), ":")) {
421   -
422   - // Get name components
423   - std::vector<std::string> name_components = gammalib::split(par.name(), ":");
424   -
425   -
426   - // Join new name
427   - std::string new_name = name_components[1];
428   - for (unsigned int i = 2; i < name_components.size(); i++) {
429   - new_name.append(":");
430   - new_name.append(name_components[i]);
431   - }
432   -
433   - // Set new name
434   - par.name(new_name);
  410 + // Get model parameter
  411 + GModelPar& par = (*spatial)[k];
  412 +
  413 + // If name contains colon then strip the prefix and the colon from
  414 + // the name
  415 + std::string name = par.name();
  416 + int found = name.find(":");
  417 + if (found != std::string::npos) {
  418 + name = name.substr(found+1, std::string::npos);
  419 + par.name(name);
435 420 }
436 421  
437 422 } // endfor looped over all parameters
438 423  
439   - // Write component
440   - if (matching_model) {
441   - // Use existing xml element
442   - spat->write(*matching_model);
443   - } else {
444   - // Create new xml element
445   - GXmlElement new_element("spatialModel");
446   - spat->write(new_element);
447   - xml.append(new_element);
  424 + // If an XML element exists then use it ...
  425 + if (matching_model != NULL) {
  426 + spatial->write(*matching_model);
448 427 }
449   -
  428 +
  429 + // ... otherwise create new XML element
  430 + else {
  431 + GXmlElement element("spatialModel");
  432 + spatial->write(element);
  433 + xml.append(element);
  434 + }
  435 +
  436 + // Delete clone
  437 + delete spatial;
  438 +
450 439 } // endfor: loop over all components
451 440  
452 441 // TODO: Remove unused xml elements
... ... @@ -457,46 +446,36 @@ void GModelSpatialComposite::write(GXmlElement&amp; xml) const
457 446  
458 447  
459 448 /***********************************************************************//**
460   - * @brief Return number of model components
461   - *
462   - * @return Number of model components
463   - ***************************************************************************/
464   -int GModelSpatialComposite::components(void) const {
465   - return m_components.size();
466   -}
467   -
468   -
469   -/***********************************************************************//**
470 449 * @brief Append spatial component
471 450 *
472   - * @param[in] model Spatial model component to append
  451 + * @param[in] component Spatial model component to append
473 452 * @param[in] name Name of spatial model (can be empty)
474 453 *
475 454 * Appends a spatial component to the composite model
476 455 ***************************************************************************/
477 456 void GModelSpatialComposite::append(const GModelSpatial& component,
478   - const std::string& name) {
479   -
480   -// Append model container
  457 + const std::string& name)
  458 +{
  459 + // Append model container
481 460 m_components.push_back(component.clone());
482 461  
483 462 // Get index of latest model
484 463 int index = m_components.size()-1;
485 464  
486 465 // Use model index if component name is empty
487   - std::string component_name = !name.empty() ? name : gammalib::str(m_components.size());
  466 + std::string component_name = !name.empty() ? name
  467 + : gammalib::str(m_components.size());
488 468  
489 469 // Check if component name is unique, throw exception if not
490 470 if (gammalib::contains(m_component_names, component_name)) {
491   - std::string msg =
492   - "Attempt to append component with name \""+component_name+"\" to composite"
493   - " spectral model container, but a component with the same name exists already.\n"
494   - "Every component in the container needs a unique name. On default the system will"
495   - "increment an integer if no component name is provided.";
  471 + std::string msg = "Attempt to append component \""+component_name+"\" "
  472 + "to composite spatial model, but a component with the "
  473 + "same name exists already. Each component needs a "
  474 + "unique name.";
496 475 throw GException::invalid_value(G_APPEND, msg);
497 476 }
498 477  
499   - // Add component name (for now simple number)
  478 + // Add component name
500 479 m_component_names.push_back(component_name);
501 480  
502 481 // Get number of spectral parameters from model
... ... @@ -515,56 +494,33 @@ void GModelSpatialComposite::append(const GModelSpatial&amp; component,
515 494 m_pars.push_back(par);
516 495  
517 496 } // endfor: loop over model parameters
518   -}
519   -
520   -std::string GModelSpatialComposite::print(const GChatter& chatter) const
521   -{
522   - // Initialise result string
523   - std::string result;
524   -
525   - // Continue only if chatter is not silent
526   - if (chatter != SILENT) {
527   -
528   - // Append header
529   - result.append("=== GModelSpatialComposite ===");
530 497  
531   - // Append information
532   - result.append("\n"+gammalib::parformat("Number of components"));
533   - result.append(gammalib::str(components()));
534   - result.append("\n"+gammalib::parformat("Number of parameters"));
535   - result.append(gammalib::str(size()));
536   -
537   - // Print parameter information
538   - for (int i = 0; i < size(); ++i) {
539   - result.append("\n"+m_pars[i]->print(chatter));
540   - }
541   -
542   - } // endif: chatter was not silent
543   -
544   - // Return result
545   - return result;
  498 + // Return
  499 + return;
546 500 }
547 501  
548 502  
549   -
550 503 /***********************************************************************//**
551 504 * @brief Returns pointer to spatial component element
552 505 *
553   - * @param[in] index Index of spatial component.
  506 + * @param[in] index Index of spatial component [0,...,components()-1].
554 507 * @return Spatial model.
555 508 *
  509 + * @exception GException::out_of_range
  510 + * Index is out of range.
  511 + *
556 512 * Returns a spatial component to the composite model
557 513 ***************************************************************************/
558 514 const GModelSpatial* GModelSpatialComposite::component(const int& index) const
559 515 {
560 516 // Check if index is in validity range
561   - if (index >= m_components.size() || index < 0) {
562   - throw GException::out_of_range(G_COMPONENT_INDEX, "Component Index", index, m_components.size(),"");
  517 + if (index < 0 || index >= m_components.size()) {
  518 + throw GException::out_of_range(G_COMPONENT_INDEX, "Component Index",
  519 + index, m_components.size());
563 520 }
564 521  
565   - // Return spectral component
  522 + // Return spatial component
566 523 return m_components[index];
567   -
568 524 }
569 525  
570 526  
... ... @@ -574,14 +530,17 @@ const GModelSpatial* GModelSpatialComposite::component(const int&amp; index) const
574 530 * @param[in] name Name of spatial component.
575 531 * @return Spatial model.
576 532 *
  533 + * @exception GException::invalid_argument
  534 + * Spatial model not found.
  535 + *
577 536 * Returns a spatial component of the composite model
578 537 ***************************************************************************/
579 538 const GModelSpatial* GModelSpatialComposite::component(const std::string& name) const
580 539 {
581 540 // Check if model name is found
582 541 int index = -1;
583   - for(int i = 0; i < m_component_names.size(); ++i) {
584   - if (m_component_names[i] == name) {
  542 + for (int i = 0; i < m_component_names.size(); ++i) {
  543 + if (m_component_names[i] == name) {
585 544 index = i;
586 545 break;
587 546 }
... ... @@ -589,14 +548,51 @@ const GModelSpatial* GModelSpatialComposite::component(const std::string&amp; name)
589 548  
590 549 // Check if component name was found
591 550 if (index == -1) {
592   - throw GException::model_not_found(G_COMPONENT_NAME, name,"");
  551 + std::string msg = "Model component \""+name+"\" not found in composite "
  552 + "spatial model.";
  553 + throw GException::invalid_argument(G_COMPONENT_NAME, msg);
593 554 }
594 555  
595   - // Return spectral component
  556 + // Return spatial component
596 557 return m_components[index];
  558 +}
  559 +
  560 +
  561 +/***********************************************************************//**
  562 + * @brief Print composite spatial model information
  563 + *
  564 + * @param[in] chatter Chattiness.
  565 + * @return String containing composite spatial model information.
  566 + ***************************************************************************/
  567 +std::string GModelSpatialComposite::print(const GChatter& chatter) const
  568 +{
  569 + // Initialise result string
  570 + std::string result;
597 571  
  572 + // Continue only if chatter is not silent
  573 + if (chatter != SILENT) {
  574 +
  575 + // Append header
  576 + result.append("=== GModelSpatialComposite ===");
  577 +
  578 + // Append information
  579 + result.append("\n"+gammalib::parformat("Number of components"));
  580 + result.append(gammalib::str(components()));
  581 + result.append("\n"+gammalib::parformat("Number of parameters"));
  582 + result.append(gammalib::str(size()));
  583 +
  584 + // Print parameter information
  585 + for (int i = 0; i < size(); ++i) {
  586 + result.append("\n"+m_pars[i]->print(chatter));
  587 + }
  588 +
  589 + } // endif: chatter was not silent
  590 +
  591 + // Return result
  592 + return result;
598 593 }
599 594  
  595 +
600 596 /*==========================================================================
601 597 = =
602 598 = Private methods =
... ... @@ -628,26 +624,29 @@ void GModelSpatialComposite::init_members(void)
628 624 void GModelSpatialComposite::copy_members(const GModelSpatialComposite& model)
629 625 {
630 626 // Copy members
631   - m_type = model.m_type;
  627 + m_type = model.m_type;
632 628 m_component_names = model.m_component_names;
  629 +
  630 + // Initialise components
633 631 m_components.clear();
634 632 m_pars.clear();
635   - for (unsigned int idx_comp = 0; idx_comp < model.m_components.size(); idx_comp ++) {
636   - m_components.push_back(model.m_components[idx_comp]->clone());
  633 +
  634 + // Copy components
  635 + for (int i = 0; i < model.m_components.size(); ++i) {
  636 +
  637 + // Clone component
  638 + m_components.push_back(model.m_components[i]->clone());
637 639  
638 640 // Get number of parameters to append
639   - int npars = m_components[idx_comp]->size();
  641 + int npars = m_components[i]->size();
640 642  
641   - // Append parameters
642   - for (unsigned int idx_par = 0; idx_par < npars; idx_par++) {
643   - // Get model parameter
644   - GModelPar& par = (*m_components[idx_comp])[idx_par];
645   -
646   - // Append to paramaters
  643 + // Append parameter references
  644 + for (int ipar = 0; ipar < npars; ++ipar) {
  645 + GModelPar& par = (*m_components[i])[ipar];
647 646 m_pars.push_back(&par);
648 647 }
649   - }
650   -
  648 +
  649 + } // endfor: looped over all components
651 650  
652 651 // Return
653 652 return;
... ... @@ -660,13 +659,14 @@ void GModelSpatialComposite::copy_members(const GModelSpatialComposite&amp; model)
660 659 void GModelSpatialComposite::free_members(void)
661 660 {
662 661 // Free model components
663   - for (unsigned int i = 0; i < m_components.size(); i ++) {
  662 + for (int i = 0; i < m_components.size(); ++i) {
664 663  
665 664 // Delete component i
666 665 delete m_components[i];
667 666  
668 667 // Signal free pointer
669 668 m_components[i] = NULL;
  669 +
670 670 }
671 671  
672 672 // Return
... ...
src/model/GModelSpatialEllipticalDisk.cpp
  View file @8947020
... ... @@ -599,7 +599,7 @@ void GModelSpatialEllipticalDisk::write(GXmlElement&amp; xml) const
599 599 /***********************************************************************//**
600 600 * @brief Print information
601 601 *
602   - * @param[in] chatter Chattiness (defaults to NORMAL).
  602 + * @param[in] chatter Chattiness.
603 603 * @return String containing model information.
604 604 ***************************************************************************/
605 605 std::string GModelSpatialEllipticalDisk::print(const GChatter& chatter) const
... ...
src/model/GModelSpatialPointSource.cpp
  View file @8947020
... ... @@ -501,6 +501,7 @@ void GModelSpatialPointSource::write(GXmlElement&amp; xml) const
501 501 /***********************************************************************//**
502 502 * @brief Print point source information
503 503 *
  504 + * @param[in] chatter Chattiness.
504 505 * @return String containing point source information.
505 506 ***************************************************************************/
506 507 std::string GModelSpatialPointSource::print(const GChatter& chatter) const
... ...
src/obs/GResponse.cpp
  View file @8947020
... ... @@ -43,8 +43,6 @@
43 43 #include "GEbounds.hpp" // will become obsolete
44 44 #include "GObservation.hpp"
45 45 #include "GModelSky.hpp"
46   -#include "GModelSpatialPointSource.hpp"
47   -#include "GModelSpatialComposite.hpp"
48 46  
49 47 /* __ Method name definitions ____________________________________________ */
50 48 #define G_IRF_RADIAL "GResponse::irf_radial(GEvent&, GSource&,"\
... ...
src/support/GTools.cpp
  View file @8947020
... ... @@ -1074,7 +1074,7 @@ bool gammalib::contains(const std::string&amp; str, const std::string&amp; substring)
1074 1074 bool result = false;
1075 1075  
1076 1076 // checks if substring is in str
1077   - if (str.find(substring) != std::string::npos){
  1077 + if (str.find(substring) != std::string::npos) {
1078 1078 result = true;
1079 1079 }
1080 1080  
... ... @@ -1092,15 +1092,18 @@ bool gammalib::contains(const std::string&amp; str, const std::string&amp; substring)
1092 1092 *
1093 1093 * Checks if a string is contained in a vector of strings
1094 1094 ***************************************************************************/
1095   -bool gammalib::contains(const std::vector<std::string> strings, const std::string& string)
  1095 +bool gammalib::contains(const std::vector<std::string> strings,
  1096 + const std::string& string)
1096 1097 {
1097 1098 // Compute result
1098   - bool result = std::find(strings.begin(), strings.end(), string) != strings.end();
  1099 + bool result = std::find(strings.begin(), strings.end(), string) !=
  1100 + strings.end();
1099 1101  
1100 1102 // Return result
1101 1103 return result;
1102 1104 }
1103 1105  
  1106 +
1104 1107 /***********************************************************************//**
1105 1108 * @brief Emits warning
1106 1109 *
... ...
test/test_GModel.cpp
  View file @8947020
... ... @@ -148,7 +148,7 @@ void TestGModel::set(void)
148 148 append(static_cast<pfunction>(&TestGModel::test_filefct),
149 149 "Test GModelSpectralFunc");
150 150 append(static_cast<pfunction>(&TestGModel::test_spectral_composite),
151   - "Test GModelSpectralComposite");
  151 + "Test GModelSpectralComposite");
152 152 append(static_cast<pfunction>(&TestGModel::test_spectral_model),
153 153 "Test spectral model XML I/O");
154 154  
... ... @@ -1410,51 +1410,42 @@ void TestGModel::test_elliptical_gauss(void)
1410 1410 ***************************************************************************/
1411 1411 void TestGModel::test_spatial_composite(void)
1412 1412 {
1413   - // Test void constructor
1414   - GModelSpatialComposite model1;
1415   - test_value(model1.type(), "Composite", "Check void model type");
1416   -
1417   - GSkyDir dir = GSkyDir();
1418   - dir.radec_deg(83.6331, 22.01);
1419   -
1420   - // Test append method
1421   - model1.append(GModelSpatialPointSource(dir));
1422   - model1.append(GModelSpatialRadialGauss(dir, 0.2));
1423   - test_value(model1.components(), 2);
1424   - test_value(model1.size(), 5);
1425   -
1426   - // Test XML constructor
1427   - GXml xml(m_xml_model_spatial_composite);
1428   - GXmlElement* element = xml.element(0)->element(0)->element("spatialModel", 0);
1429   -
1430   - GModelSpatialComposite model2(*element);
1431   - test_value(model2.size(), 5);
1432   - test_value(model2.type(), "Composite", "Check model type");
  1413 + // Test void constructor
  1414 + GModelSpatialComposite model1;
  1415 + test_value(model1.type(), "Composite", "Check void model type");
  1416 +
  1417 + // Set sky direction
  1418 + GSkyDir dir = GSkyDir();
  1419 + dir.radec_deg(83.6331, 22.01);
  1420 +
  1421 + // Test append method
  1422 + model1.append(GModelSpatialPointSource(dir));
  1423 + model1.append(GModelSpatialRadialGauss(dir, 0.2));
  1424 + test_value(model1.components(), 2);
  1425 + test_value(model1.size(), 5);
  1426 +
  1427 + // Test XML constructor
  1428 + GXml xml(m_xml_model_spatial_composite);
  1429 + GXmlElement* element = xml.element(0)->element(0)->element("spatialModel", 0);
  1430 + GModelSpatialComposite model2(*element);
  1431 + test_value(model2.size(), 5);
  1432 + test_value(model2.type(), "Composite", "Check model type");
1433 1433  
1434   - // Test access of individual parameters
1435   - test_value(model2["2:RA"].value(), 83.6331);
1436   - test_value(model2["2:DEC"].value(), 22.0145);
1437   - test_value(model2["2:Sigma"].value(), 0.2);
  1434 + // Test access of individual parameters
  1435 + test_value(model2["2:RA"].value(), 83.6331);
  1436 + test_value(model2["2:DEC"].value(), 22.0145);
  1437 + test_value(model2["2:Sigma"].value(), 0.2);
  1438 + test_value(model2["PointSource:RA"].value(), 83.6331);
  1439 + test_value(model2["2:RA"].value(), 83.6331);
  1440 + model2["2:RA"].value(83.1331);
  1441 + test_value(model2["2:RA"].value(), 83.1331);
  1442 + model2["2:DEC"].value(22.51);
  1443 + test_value(model2["2:DEC"].value(), 22.51);
  1444 + model2["2:Sigma"].value(0.6);
  1445 + test_value(model2["2:Sigma"].value(), 0.6);
1438 1446  
1439   - test_value(model2["PointSource:RA"].value(), 83.6331);
1440   -
1441   - // Test prefactor method
1442   -
1443   - test_value(model2["2:RA"].value(), 83.6331);
1444   - model2["2:RA"].value(83.1331);
1445   - test_value(model2["2:RA"].value(), 83.1331);
1446   -
1447   - // Test index method
1448   - model2["2:DEC"].value(22.51);
1449   - test_value(model2["2:DEC"].value(), 22.51);
1450   -
1451   - // Test pivot method
1452   - model2["2:Sigma"].value(0.6);
1453   - test_value(model2["2:Sigma"].value(), 0.6);
1454   -
1455   -
1456   - // Exit test
1457   - return;
  1447 + // Exit test
  1448 + return;
1458 1449 }
1459 1450  
1460 1451  
... ...