GammaLib  2.0.0
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
GCTAModelSpatialGaussSpectrum.cpp
Go to the documentation of this file.
1 /***************************************************************************
2  * GCTAModelSpatialGaussSpectrum.cpp - Spatial energy dependent Gaussian *
3  * ----------------------------------------------------------------------- *
4  * copyright (C) 2019 by Juergen Knoedlseder *
5  * ----------------------------------------------------------------------- *
6  * *
7  * This program is free software: you can redistribute it and/or modify *
8  * it under the terms of the GNU General Public License as published by *
9  * the Free Software Foundation, either version 3 of the License, or *
10  * (at your option) any later version. *
11  * *
12  * This program is distributed in the hope that it will be useful, *
13  * but WITHOUT ANY WARRANTY; without even the implied warranty of *
14  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
15  * GNU General Public License for more details. *
16  * *
17  * You should have received a copy of the GNU General Public License *
18  * along with this program. If not, see <http://www.gnu.org/licenses/>. *
19  * *
20  ***************************************************************************/
21 /**
22  * @file GCTAModelSpatialGaussSpectrum.cpp
23  * @brief Spatial energy dependent Gaussian interface implementation
24  * @author Juergen Knoedlseder
25  */
26 
27 /* __ Includes ___________________________________________________________ */
28 #ifdef HAVE_CONFIG_H
29 #include <config.h>
30 #endif
31 #include <cmath>
32 #include "GTools.hpp"
33 #include "GMath.hpp"
34 #include "GCTAInstDir.hpp"
35 #include "GEnergy.hpp"
36 #include "GTime.hpp"
40 #include "GModelSpectralConst.hpp"
41 
42 /* __ Constants __________________________________________________________ */
43 
44 /* __ Globals ____________________________________________________________ */
46 const GCTAModelSpatialRegistry g_cta_spatial_gauss_spec_registry(&g_cta_spatial_gauss_spec_seed);
47 
48 /* __ Method name definitions ____________________________________________ */
49 #define G_READ "GCTAModelSpatialGaussSpectrum::read(GXmlElement&)"
50 #define G_WRITE "GCTAModelSpatialGaussSpectrum::write(GXmlElement&)"
51 
52 /* __ Macros _____________________________________________________________ */
53 
54 /* __ Coding definitions _________________________________________________ */
55 
56 /* __ Debug definitions __________________________________________________ */
57 
58 
59 /*==========================================================================
60  = =
61  = Constructors/destructors =
62  = =
63  ==========================================================================*/
64 
65 /***********************************************************************//**
66  * @brief Void constructor
67  ***************************************************************************/
70 {
71  // Initialise members
72  init_members();
73 
74  // Return
75  return;
76 }
77 
78 
79 /***********************************************************************//**
80  * @brief Energy-independent Gaussian constructor
81  *
82  * @param[in] sigma Gaussian width (degrees\f$^2\f$).
83  ***************************************************************************/
86 {
87  // Initialise members
88  init_members();
89 
90  // Assign sigma
91  this->sigma(sigma);
92 
93  // Return
94  return;
95 }
96 
97 
98 /***********************************************************************//**
99  * @brief Gaussian spectrum constructor
100  *
101  * @param[in] sigma Spectrum defining the energy-deependent Gaussian width
102  * (degrees\f$^2\f$).
103  ***************************************************************************/
106 {
107  // Initialise members
108  init_members();
109 
110  // Assign sigma spectrum
111  this->sigma(sigma);
112 
113  // Return
114  return;
115 }
116 
117 
118 /***********************************************************************//**
119  * @brief XML constructor
120  *
121  * @param[in] xml XML element.
122  *
123  * Creates instance of a spatial energy-dependent Gaussian model by
124  * extracting information from an XML element. See
125  * GCTAModelSpatialGaussSpectrum::read() for more information about the
126  * expected structure of the XML element.
127  ***************************************************************************/
130 {
131  // Initialise members
132  init_members();
133 
134  // Read information from XML element
135  read(xml);
136 
137  // Return
138  return;
139 }
140 
141 
142 /***********************************************************************//**
143  * @brief Copy constructor
144  *
145  * @param[in] model Energy-dependent Gaussian model.
146  ***************************************************************************/
148  GCTAModelSpatial(model)
149 {
150  // Initialise members
151  init_members();
152 
153  // Copy members
154  copy_members(model);
155 
156  // Return
157  return;
158 }
159 
160 
161 /***********************************************************************//**
162  * @brief Destructor
163  ***************************************************************************/
165 {
166  // Free members
167  free_members();
168 
169  // Return
170  return;
171 }
172 
173 
174 /*==========================================================================
175  = =
176  = Operators =
177  = =
178  ==========================================================================*/
179 
180 /***********************************************************************//**
181  * @brief Assignment operator
182  *
183  * @param[in] model Energy-dependent Gaussian model.
184  ***************************************************************************/
186 {
187  // Execute only if object is not identical
188  if (this != &model) {
189 
190  // Copy base class members
191  this->GCTAModelSpatial::operator=(model);
192 
193  // Free members
194  free_members();
195 
196  // Initialise members
197  init_members();
198 
199  // Copy members
200  copy_members(model);
201 
202  } // endif: object was not identical
203 
204  // Return
205  return *this;
206 }
207 
208 
209 /*==========================================================================
210  = =
211  = Public methods =
212  = =
213  ==========================================================================*/
214 
215 /***********************************************************************//**
216  * @brief Clear instance
217  ***************************************************************************/
219 {
220  // Free class members (base and derived classes, derived class first)
221  free_members();
223 
224  // Initialise members
226  init_members();
227 
228  // Return
229  return;
230 }
231 
232 
233 /***********************************************************************//**
234  * @brief Clone instance
235  ***************************************************************************/
237 {
238  return new GCTAModelSpatialGaussSpectrum(*this);
239 }
240 
241 
242 /***********************************************************************//**
243  * @brief Evaluate function
244  *
245  * @param[in] dir Event direction.
246  * @param[in] energy Event energy.
247  * @param[in] time Event time.
248  * @param[in] gradients Compute gradients?
249  * @return Function value
250  *
251  * Evaluates the energy-dependent Gaussian model for a given event direction.
252  * The Gaussian model is defined by
253  *
254  * \f[f(\theta,E) = \exp \left(-\frac{1}{2}
255  * \left( \frac{\theta^2}{\sigma(E)} \right)^2 \right)\f]
256  *
257  * where
258  * \f$\theta\f$ is the offset angle (in degrees), and
259  * \f$\sigma(E)\f$ is the energy-dependent Gaussian width (in degrees\f$^2\f$).
260  *
261  * If the @p gradients flag is true the method will also compute the partial
262  * derivatives of the parameters. The partial derivative of the Gaussian width
263  * is given by
264  *
265  * \f[\frac{df}{d\sigma_v} = f(\theta) \frac{\theta^4}{\sigma^3} \sigma_s\f]
266  *
267  * where
268  * \f$\sigma_v\f$ is the value part,
269  * \f$\sigma_s\f$ is the scaling part, and
270  * \f$\sigma = \sigma_v \sigma_s\f$.
271  *
272  * Note that this method implements a function which is unity for
273  * \f$\theta=0\f$.
274  ***************************************************************************/
276  const GEnergy& energy,
277  const GTime& time,
278  const bool& gradients) const
279 {
280  // Initialise value
281  double value = 0.0;
282 
283  // Continue only if a sigma spectrum exists
284  if (m_sigma != NULL) {
285 
286  // Compute offset angle in degrees
287  double offset = dir.theta() * gammalib::rad2deg;
288 
289  // Compute energy dependent sigma value in degrees^2
290  double sigma = m_sigma->eval(energy, time, gradients);
291 
292  // Compute Gaussian value
293  double arg = offset * offset / sigma;
294  double arg2 = arg * arg;
295  value = std::exp(-0.5 * arg2);
296 
297  // Optionally compute partial derivatives
298  if (gradients) {
299 
300  // Compute partial derivatives of the sigma parameter
301  double g_sigma = value * arg2 / sigma;
302 
303  // Loop over model parameters
304  for (int i = 0; i < m_sigma->size(); ++i) {
305 
306  // Get reference to model parameter
307  GModelPar& par = m_sigma->operator[](i);
308 
309  // Scale parameter gradient
310  par.gradient(par.gradient()*g_sigma);
311 
312  } // endfor: loop over model parameters
313 
314  } // endif: computed partial derivatives
315 
316  // Compile option: Check for NaN/Inf
317  #if defined(G_NAN_CHECK)
318  if (gammalib::is_notanumber(value) || gammalib::is_infinite(value)) {
319  std::cout << "*** ERROR: GCTAModelSpatialGaussSpectrum::eval";
320  std::cout << "(offset=" << offset << "): NaN/Inf encountered";
321  std::cout << " (value=" << value;
322  std::cout << ", energy=" << energy;
323  std::cout << ", sigma=" << sigma;
324  std::cout << ")" << std::endl;
325  }
326  #endif
327 
328  } // endif: sigma spectrum existed
329 
330  // Return value
331  return value;
332 }
333 
334 
335 /***********************************************************************//**
336  * @brief Read model from XML element
337  *
338  * @param[in] xml XML element.
339  *
340  * Read the energy-dependent Gaussian spatial model information from an XML
341  * element. The XML element needs to be of the following format:
342  *
343  * <spatialModel type="EnergyDependentGaussian">
344  * <sigma type="...">
345  * ...
346  * </sigma>
347  * </spatialModel>
348  *
349  * where any spectral model type can be provided for the @a sigma tag.
350  ***************************************************************************/
352 {
353  // Clear model
354  clear();
355 
356  // Get pointers on sigma model components
357  const GXmlElement* sigma = xml.element("sigma", 0);
358 
359  // Get sigma model
360  GModelSpectralRegistry registry;
361  m_sigma = registry.alloc(*sigma);
362 
363  // Set pointers
364  set_pointers();
365 
366  // Return
367  return;
368 }
369 
370 
371 /***********************************************************************//**
372  * @brief Write model into XML element
373  *
374  * @param[in] xml XML element.
375  *
376  * @exception GException::invalid_value
377  * Spatial model is not of valid type.
378  *
379  * Write the energy-dependent Gaussian spatial model information into an XML
380  * element. The XML element will be of the following format:
381  *
382  * <spatialModel type="EnergyDependentGaussian">
383  * <sigma type="...">
384  * ...
385  * </sigma>
386  * </spatialModel>
387  *
388  * where the spectral sigma model will be written into the @a sigma tag.
389  ***************************************************************************/
391 {
392  // Set model type
393  if (xml.attribute("type") == "") {
394  xml.attribute("type", type());
395  }
396 
397  // Verify model type
398  if (xml.attribute("type") != type()) {
399  std::string msg = "Spatial model \""+xml.attribute("type")+
400  "\" is not of type \""+type()+"\".";
402  }
403 
404  // Write sigma spectrum if it exists
405  if (m_sigma != NULL) {
406 
407  // Create new sigma spectrum node
408  xml.append(GXmlElement("sigma"));
409 
410  // Get new sigma spectrum node
411  GXmlElement* sigma = xml.element("sigma", xml.elements("sigma")-1);
412 
413  // Write sigma spectrum
414  m_sigma->write(*sigma);
415 
416  } // endif: sigma spectrum was not NULL
417 
418  // Return
419  return;
420 }
421 
422 
423 /***********************************************************************//**
424  * @brief Set sigma spectrum from value
425  *
426  * @param[in] sigma Constant sigma value.
427  ***************************************************************************/
428 void GCTAModelSpatialGaussSpectrum::sigma(const double& sigma)
429 {
430  // Clear sigma spectrum
431  clear();
432 
433  // Set constant spectrum
434  m_sigma = new GModelSpectralConst(sigma);
435 
436  // Set parameter pointers
437  set_pointers();
438 
439  // Return
440  return;
441 }
442 
443 
444 /***********************************************************************//**
445  * @brief Set sigma spectrum
446  *
447  * @param[in] sigma Sigma spectrum.
448  ***************************************************************************/
450 {
451  // Clear sigma spectrum
452  clear();
453 
454  // Set spectrum
455  m_sigma = sigma.clone();
456 
457  // Set parameter pointers
458  set_pointers();
459 
460  // Return
461  return;
462 }
463 
464 
465 /***********************************************************************//**
466  * @brief Print point source information
467  *
468  * @param[in] chatter Chattiness.
469  * @return String containing point source information.
470  ***************************************************************************/
471 std::string GCTAModelSpatialGaussSpectrum::print(const GChatter& chatter) const
472 {
473  // Initialise result string
474  std::string result;
475 
476  // Continue only if chatter is not silent
477  if (chatter != SILENT) {
478 
479  // Append header
480  result.append("=== GCTAModelSpatialGaussSpectrum ===");
481 
482  // Determine number of sigma parameters
483  int n_sigma = (m_sigma != NULL) ? m_sigma->size() : 0;
484 
485  // Append sigma spectrum type
486  result.append("\n"+gammalib::parformat("Sigma spectrum"));
487  if (n_sigma > 0) {
488  result.append("\""+m_sigma->type()+"\"");
489  }
490 
491  // Append sigma spectrum parameters
492  for (int i = 0; i < n_sigma; ++i) {
493  result.append("\n"+(m_sigma)[i].print());
494  }
495 
496  } // endif: chatter was not silent
497 
498  // Return result
499  return result;
500 }
501 
502 
503 /*==========================================================================
504  = =
505  = Private methods =
506  = =
507  ==========================================================================*/
508 
509 /***********************************************************************//**
510  * @brief Initialise class members
511  ***************************************************************************/
513 {
514  // Initialise sigma spectrum pointer
515  m_sigma = NULL;
516 
517  // Clear parameter pointer(s)
518  m_pars.clear();
519 
520  // Return
521  return;
522 }
523 
524 
525 /***********************************************************************//**
526  * @brief Copy class members
527  *
528  * @param[in] model Energy-dependent Gaussian model.
529  ***************************************************************************/
531 {
532  // Clone sigma spectrum
533  if (model.m_sigma != NULL) {
534  m_sigma = model.m_sigma->clone();
535  }
536 
537  // Set parameter pointers
538  set_pointers();
539 
540  // Return
541  return;
542 }
543 
544 
545 /***********************************************************************//**
546  * @brief Delete class members
547  ***************************************************************************/
549 {
550  // Delete sigma spectrum
551  if (m_sigma != NULL) delete m_sigma;
552 
553  // Signal free pointer
554  m_sigma = NULL;
555 
556  // Return
557  return;
558 }
559 
560 
561 /***********************************************************************//**
562  * @brief Set pointers
563  *
564  * Set pointers to all model parameters. The pointers are stored in a vector
565  * that is member of the GModelData base class.
566  ***************************************************************************/
568 {
569  // Clear parameters
570  m_pars.clear();
571 
572  // Determine number of sigma parameters
573  int n_sigma = (m_sigma != NULL) ? m_sigma->size() : 0;
574 
575  // Gather sigma parameters if there are some
576  if (n_sigma > 0) {
577  for (int i = 0; i < n_sigma; ++i) {
578  m_pars.push_back(&((*m_sigma)[i]));
579  }
580  }
581 
582  // Return
583  return;
584 }
virtual void write(GXmlElement &xml) const =0
virtual ~GCTAModelSpatialGaussSpectrum(void)
Destructor.
Abstract spatial model class.
virtual void clear(void)
Clear instance.
Energy value class definition.
Abstract spectral model base class.
void free_members(void)
Delete class members.
double gradient(void) const
Return parameter gradient.
int size(void) const
Return number of parameters.
void init_members(void)
Initialise class members.
virtual double eval(const GEnergy &srcEng, const GTime &srcTime=GTime(), const bool &gradients=false) const =0
virtual GCTAModelSpatial & operator=(const GCTAModelSpatial &model)
Assignment operator.
XML element node class.
Definition: GXmlElement.hpp:48
Spectral model registry class definition.
virtual GCTAModelSpatialGaussSpectrum * clone(void) const
Clone instance.
virtual int elements(void) const
Return number of GXMLElement children of node.
Definition: GXmlNode.cpp:586
Time class.
Definition: GTime.hpp:55
Gammalib tools definition.
virtual std::string type(void) const =0
Interface definition for the spatial model registry class.
virtual void write(GXmlElement &xml) const
Write model into XML element.
bool is_notanumber(const double &x)
Signal if argument is not a number.
Definition: GTools.hpp:201
bool is_infinite(const double &x)
Signal if argument is infinite.
Definition: GTools.hpp:184
Model parameter class.
Definition: GModelPar.hpp:87
Spatial energy dependent Gaussian interface definition.
const GXmlAttribute * attribute(const int &index) const
Return attribute.
double theta(void) const
Return offset angle (in radians)
void free_members(void)
Delete class members.
virtual GCTAModelSpatialGaussSpectrum & operator=(const GCTAModelSpatialGaussSpectrum &model)
Assignment operator.
Constant spectral model class interface definition.
virtual std::string type(void) const
Return model type.
void init_members(void)
Initialise class members.
const GModelSpectral * sigma(void) const
Return pointer to sigma spectrum.
CTA instrument direction class interface definition.
GModelSpectral * alloc(const GXmlElement &xml) const
Allocate spectral model that is found in XML element.
GChatter
Definition: GTypemaps.hpp:33
Spatial energy dependent Gaussian model class.
Interface definition for the spectral model registry class.
void copy_members(const GCTAModelSpatialGaussSpectrum &model)
Copy class members.
const GCTAModelSpatialGaussSpectrum g_cta_spatial_gauss_spec_seed
virtual GModelSpectral * clone(void) const =0
Clones object.
virtual double eval(const GCTAInstDir &dir, const GEnergy &energy, const GTime &time, const bool &gradients=false) const
Evaluate function.
GCTAModelSpatialGaussSpectrum(void)
Void constructor.
virtual GXmlElement * element(const int &index)
Return pointer to GXMLElement child.
Definition: GXmlNode.cpp:640
virtual std::string print(const GChatter &chatter=NORMAL) const
Print point source information.
GVector exp(const GVector &vector)
Computes exponential of vector elements.
Definition: GVector.cpp:1232
virtual void read(const GXmlElement &xml)
Read model from XML element.
CTA instrument direction class.
Definition: GCTAInstDir.hpp:63
const double rad2deg
Definition: GMath.hpp:44
virtual GXmlNode * append(const GXmlNode &node)
Append XML child node.
Definition: GXmlNode.cpp:287
std::string parformat(const std::string &s, const int &indent=0)
Convert string in parameter format.
Definition: GTools.cpp:1143
std::vector< GModelPar * > m_pars
Parameter pointers.
Time class interface definition.
Spatial model registry class definition.
Constant spectral model class.
Mathematical function definitions.
Class that handles energies in a unit independent way.
Definition: GEnergy.hpp:48