GammaLib  1.7.0.dev
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
GModelSpatialRadialGauss.cpp
Go to the documentation of this file.
1 /***************************************************************************
2  * GModelSpatialRadialGauss.cpp - Radial Gaussian source model class *
3  * ----------------------------------------------------------------------- *
4  * copyright (C) 2011-2018 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 GModelSpatialRadialGauss.cpp
23  * @brief Radial Gaussian model class implementation
24  * @author Juergen Knoedlseder
25  */
26 
27 /* __ Includes ___________________________________________________________ */
28 #ifdef HAVE_CONFIG_H
29 #include <config.h>
30 #endif
31 #include "GException.hpp"
32 #include "GTools.hpp"
33 #include "GMath.hpp"
36 
37 /* __ Constants __________________________________________________________ */
38 
39 /* __ Globals ____________________________________________________________ */
41 const GModelSpatialRegistry g_radial_gauss_registry(&g_radial_gauss_seed);
42 #if defined(G_LEGACY_XML_FORMAT)
43 const GModelSpatialRadialGauss g_radial_gauss_legacy_seed(true, "GaussFunction");
44 const GModelSpatialRegistry g_radial_gauss_legacy_registry(&g_radial_gauss_legacy_seed);
45 #endif
46 
47 /* __ Method name definitions ____________________________________________ */
48 #define G_READ "GModelSpatialRadialGauss::read(GXmlElement&)"
49 #define G_WRITE "GModelSpatialRadialGauss::write(GXmlElement&)"
50 
51 /* __ Macros _____________________________________________________________ */
52 
53 /* __ Coding definitions _________________________________________________ */
54 
55 /* __ Debug definitions __________________________________________________ */
56 
57 
58 /*==========================================================================
59  = =
60  = Constructors/destructors =
61  = =
62  ==========================================================================*/
63 
64 /***********************************************************************//**
65  * @brief Void constructor
66  *
67  * Constructs empty radial Gaussian model.
68  ***************************************************************************/
70 {
71  // Initialise members
72  init_members();
73 
74  // Return
75  return;
76 }
77 
78 
79 /***********************************************************************//**
80  * @brief Model type constructor
81  *
82  * @param[in] dummy Dummy flag.
83  * @param[in] type Model type.
84  *
85  * Constructs empty radial Gaussian model by specifying a model @p type.
86  ***************************************************************************/
88  const std::string& type) :
90 {
91  // Initialise members
92  init_members();
93 
94  // Set model type
95  m_type = type;
96 
97  // Return
98  return;
99 }
100 
101 
102 /***********************************************************************//**
103  * @brief Constructor
104  *
105  * @param[in] dir Sky position of Gaussian.
106  * @param[in] sigma Width of Gaussian (in degrees).
107  *
108  * Constructs a Gaussian spatial model using a sky direction (@p dir) and
109  * a Gaussian width parameter @p sigma in degrees.
110  ***************************************************************************/
112  const double& sigma) :
114 {
115  // Initialise members
116  init_members();
117 
118  // Assign direction and sigma
119  this->dir(dir);
120  this->sigma(sigma);
121 
122  // Return
123  return;
124 }
125 
126 
127 /***********************************************************************//**
128  * @brief XML constructor
129  *
130  * @param[in] xml XML element.
131  *
132  * Constructs a Gaussian spatial model by extracting information from an XML
133  * element. See the method read() for more information about the expected
134  * structure of the XML element.
135  ***************************************************************************/
138 {
139  // Initialise members
140  init_members();
141 
142  // Read information from XML element
143  read(xml);
144 
145  // Return
146  return;
147 }
148 
149 
150 /***********************************************************************//**
151  * @brief Copy constructor
152  *
153  * @param[in] model Radial Gaussian model.
154  ***************************************************************************/
156  GModelSpatialRadial(model)
157 {
158  // Initialise members
159  init_members();
160 
161  // Copy members
162  copy_members(model);
163 
164  // Return
165  return;
166 }
167 
168 
169 /***********************************************************************//**
170  * @brief Destructor
171  ***************************************************************************/
173 {
174  // Free members
175  free_members();
176 
177  // Return
178  return;
179 }
180 
181 
182 /*==========================================================================
183  = =
184  = Operators =
185  = =
186  ==========================================================================*/
187 
188 /***********************************************************************//**
189  * @brief Assignment operator
190  *
191  * @param[in] model Radial Gaussian model.
192  * @return Radial Gaussian model.
193  ***************************************************************************/
195 {
196  // Execute only if object is not identical
197  if (this != &model) {
198 
199  // Copy base class members
200  this->GModelSpatialRadial::operator=(model);
201 
202  // Free members
203  free_members();
204 
205  // Initialise members
206  init_members();
207 
208  // Copy members
209  copy_members(model);
210 
211  } // endif: object was not identical
212 
213  // Return
214  return *this;
215 }
216 
217 
218 /*==========================================================================
219  = =
220  = Public methods =
221  = =
222  ==========================================================================*/
223 
224 /***********************************************************************//**
225  * @brief Clear radial Gauss model
226  ***************************************************************************/
228 {
229  // Free class members (base and derived classes, derived class first)
230  free_members();
233 
234  // Initialise members
237  init_members();
238 
239  // Return
240  return;
241 }
242 
243 
244 /***********************************************************************//**
245  * @brief Clone radial Gauss model
246  ***************************************************************************/
248 {
249  // Clone radial Gauss model
250  return new GModelSpatialRadialGauss(*this);
251 }
252 
253 
254 /***********************************************************************//**
255  * @brief Evaluate function
256  *
257  * @param[in] theta Angular distance from Gaussian centre (radians).
258  * @param[in] energy Photon energy.
259  * @param[in] time Photon arrival time.
260  * @param[in] gradients Compute gradients?
261  * @return Model value.
262  *
263  * Evaluates the spatial component for a Gaussian source model using
264  *
265  * \f[
266  * S_{\rm p}(\vec{p} | E, t) =
267  * \frac{1}{2 \pi \sigma^2} \exp
268  * \left(-\frac{1}{2}\frac{\theta^2}{\sigma^2} \right)
269  * \f]
270  *
271  * where
272  * - \f$\theta\f$ is the angular separation from the centre of the model, and
273  * - \f$\sigma\f$ is the Gaussian width.
274  *
275  * The method will not compute analytical parameter gradients, even if the
276  * @p gradients argument is set to true. Radial disk parameter gradients
277  * need to be computed numerically.
278  *
279  * @todo The Gaussian function is only correct in the small angle
280  * approximation.
281  ***************************************************************************/
282 double GModelSpatialRadialGauss::eval(const double& theta,
283  const GEnergy& energy,
284  const GTime& time,
285  const bool& gradients) const
286 {
287  // Compute value
288  double sigma_rad = sigma() * gammalib::deg2rad;
289  double sigma2 = sigma_rad * sigma_rad;
290  double theta2 = theta * theta;
291  double value = std::exp(-0.5 * theta2 / sigma2) /
292  (gammalib::twopi * sigma2);
293 
294  // Compile option: Check for NaN/Inf
295  #if defined(G_NAN_CHECK)
296  if (gammalib::is_notanumber(value) || gammalib::is_infinite(value)) {
297  std::cout << "*** ERROR: GModelSpatialRadialGauss::eval";
298  std::cout << "(theta=" << theta << "): NaN/Inf encountered";
299  std::cout << " (value=" << value;
300  std::cout << ", sigma_rad=" << sigma_rad;
301  std::cout << ", sigma2=" << sigma2;
302  std::cout << ", theta2=" << theta2;
303  std::cout << ")" << std::endl;
304  }
305  #endif
306 
307  // Return value
308  return value;
309 }
310 
311 
312 /***********************************************************************//**
313  * @brief Returns MC sky direction
314  *
315  * @param[in] energy Photon energy.
316  * @param[in] time Photon arrival time.
317  * @param[in,out] ran Random number generator.
318  * @return Sky direction.
319  *
320  * Draws an arbitrary sky direction from the 2D Gaussian distribution as
321  * function of the photon @p energy and arrival @p time.
322  *
323  * @todo This method is only valid in the small angle approximation.
324  ***************************************************************************/
326  const GTime& time,
327  GRan& ran) const
328 {
329  // Simulate offset from photon arrival direction
330  double theta = sigma() * ran.chisq2();
331  double phi = 360.0 * ran.uniform();
332 
333  // Rotate sky direction by offset
334  GSkyDir sky_dir = dir();
335  sky_dir.rotate_deg(phi, theta);
336 
337  // Return sky direction
338  return sky_dir;
339 }
340 
341 
342 /***********************************************************************//**
343  * @brief Checks where model contains specified sky direction
344  *
345  * @param[in] dir Sky direction.
346  * @param[in] margin Margin to be added to sky direction (degrees)
347  * @return True if the model contains the sky direction.
348  *
349  * Signals whether a sky direction is contained in the radial gauss model.
350  ***************************************************************************/
352  const double& margin) const
353 {
354  // Compute distance to centre (radians)
355  double distance = dir.dist(this->dir());
356 
357  // Return flag
358  return (distance <= theta_max() + margin*gammalib::deg2rad);
359 }
360 
361 
362 /***********************************************************************//**
363  * @brief Return maximum model radius (in radians)
364  *
365  * Returns \f$5 \sigma\f$ as approximate edge of the Gaussian. This limit
366  * is of course arbitrary, but allows to limit the integration region for
367  * response computation.
368  ***************************************************************************/
370 {
371  // Return value
372  return (sigma() * gammalib::deg2rad * 5.0);
373 }
374 
375 
376 /***********************************************************************//**
377  * @brief Read model from XML element
378  *
379  * @param[in] xml XML element.
380  *
381  * @exception GException::model_invalid_parnum
382  * Invalid number of model parameters found in XML element.
383  * @exception GException::model_invalid_parnames
384  * Invalid model parameter names found in XML element.
385  *
386  * Reads the radial Gauss model information from an XML element. The XML
387  * element shall have either the format
388  *
389  * <spatialModel type="RadialGaussian">
390  * <parameter name="RA" scale="1.0" value="83.6331" min="-360" max="360" free="1"/>
391  * <parameter name="DEC" scale="1.0" value="22.0145" min="-90" max="90" free="1"/>
392  * <parameter name="Sigma" scale="1.0" value="0.45" min="0.01" max="10" free="1"/>
393  * </spatialModel>
394  *
395  * or
396  *
397  * <spatialModel type="RadialGaussian">
398  * <parameter name="GLON" scale="1.0" value="83.6331" min="-360" max="360" free="1"/>
399  * <parameter name="GLAT" scale="1.0" value="22.0145" min="-90" max="90" free="1"/>
400  * <parameter name="Sigma" scale="1.0" value="0.45" min="0.01" max="10" free="1"/>
401  * </spatialModel>
402  *
403  * @todo Implement a test of the sigma and sigma boundary. The sigma
404  * and sigma minimum should be >0.
405  ***************************************************************************/
407 {
408  // Determine number of parameter nodes in XML element
409  int npars = xml.elements("parameter");
410 
411  // Verify that XML element has exactly 3 parameters
412  if (xml.elements() != 3 || npars != 3) {
414  "Gaussian source model requires exactly 3 parameters.");
415  }
416 
417  // Read Gaussian location
419 
420  // Get parameters
422 
423  // Read parameters
424  m_sigma.read(*sigma);
425 
426  // Return
427  return;
428 }
429 
430 
431 /***********************************************************************//**
432  * @brief Write model into XML element
433  *
434  * @param[in] xml XML element into which model information is written.
435  *
436  * @exception GException::model_invalid_spatial
437  * Existing XML element is not of type 'GaussFunction'
438  * @exception GException::model_invalid_parnum
439  * Invalid number of model parameters found in XML element.
440  * @exception GException::model_invalid_parnames
441  * Invalid model parameter names found in XML element.
442  *
443  * Writes the radial disk model information into an XML element. The XML
444  * element will have the format
445  *
446  * <spatialModel type="RadialGaussian">
447  * <parameter name="RA" scale="1.0" value="83.6331" min="-360" max="360" free="1"/>
448  * <parameter name="DEC" scale="1.0" value="22.0145" min="-90" max="90" free="1"/>
449  * <parameter name="Sigma" scale="1.0" value="0.45" min="0.01" max="10" free="1"/>
450  * </spatialModel>
451  *
452  ***************************************************************************/
454 {
455  // Write Gaussian location
457 
458  // Get or create parameters
460 
461  // Write parameters
462  m_sigma.write(*sigma);
463 
464  // Return
465  return;
466 }
467 
468 
469 /***********************************************************************//**
470  * @brief Print Gaussian source information
471  *
472  * @param[in] chatter Chattiness (defaults to NORMAL).
473  * @return String containing model information.
474  ***************************************************************************/
475 std::string GModelSpatialRadialGauss::print(const GChatter& chatter) const
476 {
477  // Initialise result string
478  std::string result;
479 
480  // Continue only if chatter is not silent
481  if (chatter != SILENT) {
482 
483  // Append header
484  result.append("=== GModelSpatialRadialGauss ===");
485 
486  // Append parameters
487  result.append("\n"+gammalib::parformat("Number of parameters"));
488  result.append(gammalib::str(size()));
489  for (int i = 0; i < size(); ++i) {
490  result.append("\n"+m_pars[i]->print(chatter));
491  }
492 
493  } // endif: chatter was not silent
494 
495  // Return result
496  return result;
497 }
498 
499 
500 /*==========================================================================
501  = =
502  = Private methods =
503  = =
504  ==========================================================================*/
505 
506 /***********************************************************************//**
507  * @brief Initialise class members
508  *
509  * Note that this model implements no gradients, as spatial models will
510  * always require numerical gradient computations. The minimum Gaussian
511  * width is set to 1 arcsec.
512  ***************************************************************************/
514 {
515  // Initialise model type
516  m_type = "RadialGaussian";
517 
518  // Initialise Gaussian sigma
519  m_sigma.clear();
520  m_sigma.name("Sigma");
521  m_sigma.unit("deg");
522  m_sigma.value(2.778e-4); // 1 arcsec
523  m_sigma.min(2.778e-4); // 1 arcsec
524  m_sigma.free();
525  m_sigma.scale(1.0);
526  m_sigma.gradient(0.0);
527  m_sigma.has_grad(false); // Radial components never have gradients
528 
529  // Set parameter pointer(s)
530  m_pars.push_back(&m_sigma);
531 
532  // Initialise other members
533  m_region.clear();
534 
535  // Return
536  return;
537 }
538 
539 
540 /***********************************************************************//**
541  * @brief Copy class members
542  *
543  * @param[in] model Gaussian spatial model.
544  *
545  * We do not have to push back the members on the parameter stack as this
546  * should have been done by init_members() that was called before. Otherwise
547  * we would have sigma twice on the stack.
548  ***************************************************************************/
550 {
551  // Copy members
552  m_type = model.m_type;
553  m_sigma = model.m_sigma;
554  m_region = model.m_region;
555 
556  // Return
557  return;
558 }
559 
560 
561 /***********************************************************************//**
562  * @brief Delete class members
563  ***************************************************************************/
565 {
566  // Return
567  return;
568 }
569 
570 
571 /***********************************************************************//**
572  * @brief Set boundary sky region
573  ***************************************************************************/
575 {
576  // Set sky region centre to Gaussian centre
578 
579  // Set sky region radius to 5 times the Gaussian sigma
580  m_region.radius(m_sigma.value() * 5.0);
581 
582  // Return
583  return;
584 }
double chisq2(void)
Returns Chi2 deviates for 2 degrees of freedom.
Definition: GRan.cpp:370
double sigma(void) const
Return Gaussian sigma.
const std::string & name(void) const
Return parameter name.
virtual double theta_max(void) const
Return maximum model radius (in radians)
double gradient(void) const
Return parameter gradient.
int size(void) const
Return number of parameters.
void write(GXmlElement &xml) const
Set or update parameter attributes in XML element.
Definition: GModelPar.cpp:353
XML element node class.
Definition: GXmlElement.hpp:47
Random number generator class.
Definition: GRan.hpp:44
#define G_WRITE
virtual int elements(void) const
Return number of GXMLElement children of node.
Definition: GXmlNode.cpp:580
Time class.
Definition: GTime.hpp:54
Gammalib tools definition.
GModelPar m_dec
Declination (deg)
Interface definition for the spatial model registry class.
virtual void write(GXmlElement &xml) const
Write model into XML element.
double min(void) const
Return parameter minimum boundary.
void init_members(void)
Initialise class members.
GModelPar m_sigma
Gaussian width (deg)
void copy_members(const GModelSpatialRadialGauss &model)
Copy class members.
virtual GModelSpatialRadialGauss * clone(void) const
Clone radial Gauss model.
bool is_notanumber(const double &x)
Signal if argument is not a number.
Definition: GTools.hpp:184
const double & radius(void) const
Return circular region radius (in degrees)
const double & scale(void) const
Return parameter scale.
bool is_infinite(const double &x)
Signal if argument is infinite.
Definition: GTools.hpp:167
const double deg2rad
Definition: GMath.hpp:43
bool has_grad(void) const
Signal if parameter gradient is computed analytically.
virtual GSkyDir mc(const GEnergy &energy, const GTime &time, GRan &ran) const
Returns MC sky direction.
virtual std::string type(void) const
Return model type.
GModelSpatialRadialGauss(void)
Void constructor.
void free_members(void)
Delete class members.
void free(void)
Free a parameter.
std::vector< GModelPar * > m_pars
Parameter pointers.
void rotate_deg(const double &phi, const double &theta)
Rotate sky direction by zenith and azimuth angle.
Definition: GSkyDir.cpp:321
GSkyRegionCircle m_region
Bounding circle.
GXmlElement * xml_need_par(const std::string &origin, GXmlElement &xml, const std::string &name)
Return pointer to parameter with given name in XML element.
Definition: GTools.cpp:1513
virtual void read(const GXmlElement &xml)
Read model from XML element.
double uniform(void)
Returns random double precision floating value in range 0 to 1.
Definition: GRan.cpp:242
void clear(void)
Clear parameter.
Spatial model registry class definition.
virtual void read(const GXmlElement &xml)
Read model from XML element.
GChatter
Definition: GTypemaps.hpp:33
void set_region(void) const
Set boundary sky region.
virtual double eval(const double &theta, const GEnergy &energy, const GTime &time, const bool &gradients=false) const
Evaluate function.
std::string m_type
Model type.
GModelPar m_ra
Right Ascension (deg)
virtual void write(GXmlElement &xml) const
Write model into XML element.
#define G_READ
virtual GModelSpatialRadial & operator=(const GModelSpatialRadial &model)
Assignment operator.
const GSkyDir & centre(void) const
Return circular region centre.
void free_members(void)
Delete class members.
GSkyDir dir(void) const
Return position of radial spatial model.
void init_members(void)
Initialise class members.
void read(const GXmlElement &xml)
Extract parameter attributes from XML element.
Definition: GModelPar.cpp:229
double value(void) const
Return parameter value.
const std::string & unit(void) const
Return parameter unit.
Exception handler interface definition.
double dist(const GSkyDir &dir) const
Compute angular distance between sky directions in radians.
Definition: GSkyDir.hpp:265
GVector exp(const GVector &vector)
Computes exponential of vector elements.
Definition: GVector.cpp:1142
virtual GModelSpatialRadialGauss & operator=(const GModelSpatialRadialGauss &model)
Assignment operator.
Abstract radial spatial model base class.
void init_members(void)
Initialise class members.
const double twopi
Definition: GMath.hpp:36
std::string parformat(const std::string &s, const int &indent=0)
Convert string in parameter format.
Definition: GTools.cpp:1022
void free_members(void)
Delete class members.
virtual bool contains(const GSkyDir &dir, const double &margin=0.0) const
Checks where model contains specified sky direction.
Radial Gaussian model class.
Sky direction class.
Definition: GSkyDir.hpp:62
virtual void clear(void)
Clear radial Gauss model.
Radial Gaussian model class interface definition.
const GXmlElement * xml_get_par(const std::string &origin, const GXmlElement &xml, const std::string &name)
Return pointer to parameter with given name in XML element.
Definition: GTools.cpp:1562
const GModelSpatialRadialGauss g_radial_gauss_seed
Mathematical function definitions.
virtual std::string print(const GChatter &chatter=NORMAL) const
Print Gaussian source information.
void clear(void)
Clear instance.
Class that handles energies in a unit independent way.
Definition: GEnergy.hpp:48
virtual ~GModelSpatialRadialGauss(void)
Destructor.
std::string str(const unsigned short int &value)
Convert unsigned short integer value into string.
Definition: GTools.cpp:413