GammaLib  1.7.0.dev
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
GCTAModelRadialGauss.cpp
Go to the documentation of this file.
1 /***************************************************************************
2  * GCTAModelRadialGauss.cpp - Radial Gaussian CTA 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 GCTAModelRadialGauss.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 <cmath>
32 #include "GException.hpp"
33 #include "GTools.hpp"
34 #include "GMath.hpp"
35 #include "GRan.hpp"
36 #include "GIntegral.hpp"
37 #include "GCTAObservation.hpp"
38 #include "GCTAInstDir.hpp"
39 #include "GCTAModelRadialGauss.hpp"
42 
43 /* __ Constants __________________________________________________________ */
44 
45 /* __ Globals ____________________________________________________________ */
47 const GCTAModelRadialRegistry g_cta_radial_gauss_registry(&g_cta_radial_gauss_seed);
48 const GCTAModelSpatialRegistry g_cta_radial_gauss_spatial_registry(&g_cta_radial_gauss_seed);
49 
50 /* __ Method name definitions ____________________________________________ */
51 #define G_READ "GCTAModelRadialGauss::read(GXmlElement&)"
52 #define G_WRITE "GCTAModelRadialGauss::write(GXmlElement&)"
53 
54 /* __ Macros _____________________________________________________________ */
55 
56 /* __ Coding definitions _________________________________________________ */
57 //#define G_DEBUG_MC //!< Debug MC method
58 
59 /* __ Debug definitions __________________________________________________ */
60 
61 
62 /*==========================================================================
63  = =
64  = Constructors/destructors =
65  = =
66  ==========================================================================*/
67 
68 /***********************************************************************//**
69  * @brief Void constructor
70  ***************************************************************************/
72 {
73  // Initialise members
74  init_members();
75 
76  // Return
77  return;
78 }
79 
80 
81 /***********************************************************************//**
82  * @brief Constructor
83  *
84  * @param[in] sigma Gaussian width (degrees\f$^2\f$).
85  ***************************************************************************/
87 {
88  // Initialise members
89  init_members();
90 
91  // Assign sigma
92  this->sigma(sigma);
93 
94  // Return
95  return;
96 }
97 
98 
99 /***********************************************************************//**
100  * @brief XML constructor
101  *
102  * @param[in] xml XML element.
103  *
104  * Creates instance of a radial Gaussian model by extracting information
105  * from an XML element. See GCTAModelRadialGauss::read() for more information
106  * about the expected structure of the XML element.
107  ***************************************************************************/
109 {
110  // Initialise members
111  init_members();
112 
113  // Read information from XML element
114  read(xml);
115 
116  // Return
117  return;
118 }
119 
120 
121 /***********************************************************************//**
122  * @brief Copy constructor
123  *
124  * @param[in] model Radial Gaussian model.
125  ***************************************************************************/
127  GCTAModelRadial(model)
128 {
129  // Initialise members
130  init_members();
131 
132  // Copy members
133  copy_members(model);
134 
135  // Return
136  return;
137 }
138 
139 
140 /***********************************************************************//**
141  * @brief Destructor
142  ***************************************************************************/
144 {
145  // Free members
146  free_members();
147 
148  // Return
149  return;
150 }
151 
152 
153 /*==========================================================================
154  = =
155  = Operators =
156  = =
157  ==========================================================================*/
158 
159 /***********************************************************************//**
160  * @brief Assignment operator
161  *
162  * @param[in] model Radial Gaussian model.
163  ***************************************************************************/
165 {
166  // Execute only if object is not identical
167  if (this != &model) {
168 
169  // Copy base class members
170  this->GCTAModelRadial::operator=(model);
171 
172  // Free members
173  free_members();
174 
175  // Initialise members
176  init_members();
177 
178  // Copy members
179  copy_members(model);
180 
181  } // endif: object was not identical
182 
183  // Return
184  return *this;
185 }
186 
187 
188 /*==========================================================================
189  = =
190  = Public methods =
191  = =
192  ==========================================================================*/
193 
194 /***********************************************************************//**
195  * @brief Clear instance
196  ***************************************************************************/
198 {
199  // Free class members (base and derived classes, derived class first)
200  free_members();
202 
203  // Initialise members
205  init_members();
206 
207  // Return
208  return;
209 }
210 
211 
212 /***********************************************************************//**
213  * @brief Clone instance
214  ***************************************************************************/
216 {
217  return new GCTAModelRadialGauss(*this);
218 }
219 
220 
221 /***********************************************************************//**
222  * @brief Evaluate function
223  *
224  * @param[in] offset Offset angle (degrees).
225  * @param[in] gradients Compute gradients?
226  * @return Function value
227  *
228  * Evaluates the Gaussian model for a given offset. The Gaussian model is
229  * defined as
230  * \f[f(\theta) = \exp \left(-\frac{1}{2}
231  * \left( \frac{\theta^2}{\sigma} \right)^2 \right)\f]
232  * where
233  * \f$\theta\f$ is the offset angle (in degrees), and
234  * \f$\sigma\f$ is the Gaussian width (in degrees\f$^2\f$).
235  *
236  * If the @p gradients flag is true the method will also compute the partial
237  * derivatives of the parameters. The partial derivative of the Gaussian width
238  * is given by
239  * \f[\frac{df}{d\sigma_v} = f(\theta) \frac{\theta^4}{\sigma^3} \sigma_s\f]
240  * where
241  * \f$\sigma_v\f$ is the value part,
242  * \f$\sigma_s\f$ is the scaling part, and
243  * \f$\sigma = \sigma_v \sigma_s\f$.
244  *
245  * Note that this method implements a function which is unity for
246  * \f$\theta=0\f$.
247  ***************************************************************************/
248 double GCTAModelRadialGauss::eval(const double& offset,
249  const bool& gradients) const
250 {
251  // Compute value
252  double arg = offset * offset / sigma();
253  double arg2 = arg * arg;
254  double value = std::exp(-0.5 * arg2);
255 
256  // Optionally compute partial derivatives
257  if (gradients) {
258 
259  // Compute partial derivatives of the sigma parameter.
260  double g_sigma = value * arg2 / sigma() * m_sigma.scale();
261 
262  // Set factor gradient
263  m_sigma.factor_gradient(g_sigma);
264 
265  } // endif: computed partial derivatives
266 
267  // Compile option: Check for NaN/Inf
268  #if defined(G_NAN_CHECK)
269  if (gammalib::is_notanumber(value) || gammalib::is_infinite(value)) {
270  std::cout << "*** ERROR: GCTAModelRadialGauss::eval";
271  std::cout << "(offset=" << offset << "): NaN/Inf encountered";
272  std::cout << " (value=" << value;
273  std::cout << ", sigma=" << sigma();
274  std::cout << ")" << std::endl;
275  }
276  #endif
277 
278  // Return value
279  return value;
280 }
281 
282 
283 /***********************************************************************//**
284  * @brief Returns MC instrument direction
285  *
286  * @param[in,out] ran Random number generator.
287  * @return CTA instrument direction.
288  *
289  * Draws an arbitrary CTA instrument position from
290  * \f[f(\theta) = \sin(\theta)
291  * \exp \left(-\frac{1}{2}\frac{\theta^4}{\sigma^2} \right)\f]
292  * where
293  * \f$\theta\f$ is the offset angle (in degrees), and
294  * \f$\sigma\f$ is the Gaussian width (in degrees\f$^2\f$),
295  * using the rejection method.
296  *
297  * @todo Method can be optimised by using a random deviate of sin instead
298  * of the uniform random deviate which leads to many unnecessary
299  * rejections.
300  ***************************************************************************/
302 {
303  // Simulate offset from photon arrival direction
304  #if defined(G_DEBUG_MC)
305  int n_samples = 0;
306  #endif
307  double sigma_max = 4.0 * std::sqrt(sigma());
308  double u_max = sin(sigma_max * gammalib::deg2rad);
309  double value = 0.0;
310  double u = 1.0;
311  double offset = 0.0;
312  do {
313  offset = ran.uniform() * sigma_max;
314  double arg = offset * offset / sigma();
315  double arg2 = arg * arg;
316  value = sin(offset * gammalib::deg2rad) * exp(-0.5 * arg2);
317  u = ran.uniform() * u_max;
318  #if defined(G_DEBUG_MC)
319  n_samples++;
320  #endif
321  } while (u > value);
322  #if defined(G_DEBUG_MC)
323  std::cout << "#=" << n_samples << " ";
324  #endif
325 
326  // Simulate azimuth angle
327  double phi = 360.0 * ran.uniform();
328 
329  // Convert from degrees to radians
330  offset *= gammalib::deg2rad;
331  phi *= gammalib::deg2rad;
332 
333  // Compute DETX and DETY coordinates
334  double detx(0.0);
335  double dety(0.0);
336  if (offset > 0.0 ) {
337  detx = offset * std::cos(phi);
338  dety = offset * std::sin(phi);
339  }
340 
341  // Set instrument direction
342  GCTAInstDir dir(detx, dety);
343 
344  // Return instrument direction
345  return dir;
346 }
347 
348 
349 /***********************************************************************//**
350  * @brief Return maximum function value for Monte Carlo simulations
351  *
352  * @param[in] obs CTA Observation.
353  * @return Maximum function value for Monte Carlo simulations.
354  ***************************************************************************/
356 {
357  // Return maximum value
358  return 1.0;
359 }
360 
361 
362 /***********************************************************************//**
363  * @brief Returns integral over radial model (in steradians)
364  *
365  * Computes
366  * \f[\Omega = 2 \pi \int_0^{\pi} \sin \theta f(\theta) d\theta\f]
367  * where
368  * \f[f(\theta) = \exp \left(-\frac{1}{2}
369  * \left( \frac{\theta^2}{\sigma} \right)^2 \right)\f]
370  * \f$\theta\f$ is the offset angle (in degrees), and
371  * \f$\sigma\f$ is the Gaussian width (in degrees\f$^2\f$).
372  *
373  * The integration is performed numerically, and the upper integration bound
374  * \f$\pi\f$
375  * is set to
376  * \f$\sqrt(10 \sigma)\f$
377  * to reduce inaccuracies in the numerical integration.
378  ***************************************************************************/
379 double GCTAModelRadialGauss::omega(void) const
380 {
381  // Allocate integrand
385 
386  // Allocate intergal
387  GIntegral integral(&integrand);
388 
389  // Set upper integration boundary
390  double offset_max = sqrt(10.0*sigma()) * gammalib::deg2rad;
391  if (offset_max > gammalib::pi) {
392  offset_max = gammalib::pi;
393  }
394 
395  // Perform numerical integration
396  double omega = integral.romberg(0.0, offset_max) * gammalib::twopi;
397 
398  // Return integral
399  return omega;
400 }
401 
402 
403 /***********************************************************************//**
404  * @brief Read model from XML element
405  *
406  * @param[in] xml XML element.
407  *
408  * @exception GException::model_invalid_parnum
409  * Invalid number of model parameters found in XML element.
410  * @exception GException::model_invalid_parvalue
411  * Non-positive parameter value found.
412  * @exception GException::model_invalid_parlimit
413  * Missing or non-positive minimum parameter boundary.
414  * @exception GException::model_invalid_parnames
415  * Invalid model parameter names found in XML element.
416  *
417  * Read the Gaussian radial model information from an XML element. The XML
418  * element is required to have one parameter named "Sigma".
419  ***************************************************************************/
421 {
422  // Verify that XML element has exactly 1 parameter
423  if (xml.elements() != 1 || xml.elements("parameter") != 1) {
425  "Radial Gaussian model requires exactly 1 parameter.");
426  }
427 
428  // Extract model parameters
429  int npar[] = {0};
430  for (int i = 0; i < 1; ++i) {
431 
432  // Get parameter element
433  const GXmlElement* par = xml.element("parameter", i);
434 
435  // Handle Sigma
436  if (par->attribute("name") == "Sigma") {
437 
438  // Read parameter
439  m_sigma.read(*par);
440 
441  // Check parameter
442  if (m_sigma.value() <= 0.0) {
444  "\"Sigma\" parameter is required to be positive.");
445  }
446  if (!m_sigma.has_min() || m_sigma.min() <= 0.0) {
448  "\"Sigma\" parameter requires positive minimum boundary.");
449  }
450 
451  // Increment parameter counter
452  npar[0]++;
453  }
454 
455  } // endfor: looped over all parameters
456 
457  // Verify that all parameters were found
458  if (npar[0] != 1) {
460  "Radial Gaussian model requires \"Sigma\" parameter.");
461  }
462 
463  // Return
464  return;
465 }
466 
467 
468 /***********************************************************************//**
469  * @brief Write model into XML element
470  *
471  * @param[in] xml XML element.
472  *
473  * @exception GException::model_invalid_spatial
474  * Existing XML element is not of type 'GaussFunction'
475  * @exception GException::model_invalid_parnum
476  * Invalid number of model parameters found in XML element.
477  * @exception GException::model_invalid_parnames
478  * Invalid model parameter names found in XML element.
479  *
480  * Write the Gaussian radial model information into an XML element. The XML
481  * element will have one parameter leaf named "Sigma".
482  ***************************************************************************/
484 {
485  // Set model type
486  if (xml.attribute("type") == "") {
487  xml.attribute("type", type());
488  }
489 
490  // Verify model type
491  if (xml.attribute("type") != type()) {
493  "Radial Gaussian model is not of type \""+type()+"\".");
494  }
495 
496  // If XML element has 0 nodes then append 1 parameter node
497  if (xml.elements() == 0) {
498  xml.append(GXmlElement("parameter name=\"Sigma\""));
499  }
500 
501  // Verify that XML element has exactly 1 parameter
502  if (xml.elements() != 1 || xml.elements("parameter") != 1) {
504  "Radial Gaussian model requires exactly 1 parameter.");
505  }
506 
507  // Set or update model parameter attributes
508  int npar[] = {0};
509  for (int i = 0; i < 1; ++i) {
510 
511  // Get parameter element
512  GXmlElement* par = xml.element("parameter", i);
513 
514  // Handle sigma
515  if (par->attribute("name") == "Sigma") {
516  m_sigma.write(*par);
517  npar[0]++;
518  }
519 
520  } // endfor: looped over all parameters
521 
522  // Check of all required parameters are present
523  if (npar[0] != 1) {
525  "Radial Gaussian model requires \"Sigma\" parameter.");
526  }
527 
528  // Return
529  return;
530 }
531 
532 
533 /***********************************************************************//**
534  * @brief Print point source information
535  *
536  * @param[in] chatter Chattiness.
537  * @return String containing point source information.
538  ***************************************************************************/
539 std::string GCTAModelRadialGauss::print(const GChatter& chatter) const
540 {
541  // Initialise result string
542  std::string result;
543 
544  // Continue only if chatter is not silent
545  if (chatter != SILENT) {
546 
547  // Append header
548  result.append("=== GCTAModelRadialGauss ===");
549 
550  // Append information
551  result.append("\n"+gammalib::parformat("Number of parameters") +
552  gammalib::str(size()));
553  for (int i = 0; i < size(); ++i) {
554  result.append("\n"+m_pars[i]->print(chatter));
555  }
556 
557  } // endif: chatter was not silent
558 
559  // Return result
560  return result;
561 }
562 
563 
564 /*==========================================================================
565  = =
566  = Private methods =
567  = =
568  ==========================================================================*/
569 
570 /***********************************************************************//**
571  * @brief Initialise class members
572  ***************************************************************************/
574 {
575  // Initialise Gaussian sigma
576  m_sigma.clear();
577  m_sigma.name("Sigma");
578  m_sigma.unit("deg2");
579  m_sigma.value(7.71728e-8); // (1 arcsec)^2
580  m_sigma.min(7.71728e-8); // (1 arcsec)^2
581  m_sigma.free();
582  m_sigma.scale(1.0);
583  m_sigma.gradient(0.0);
584  m_sigma.has_grad(true);
585 
586  // Set parameter pointer(s)
587  m_pars.clear();
588  m_pars.push_back(&m_sigma);
589 
590  // Return
591  return;
592 }
593 
594 
595 /***********************************************************************//**
596  * @brief Copy class members
597  *
598  * @param[in] model Radial Gaussian model.
599  ***************************************************************************/
601 {
602  // Copy members
603  m_sigma = model.m_sigma;
604 
605  // Set parameter pointer(s)
606  m_pars.clear();
607  m_pars.push_back(&m_sigma);
608 
609  // Return
610  return;
611 }
612 
613 
614 /***********************************************************************//**
615  * @brief Delete class members
616  ***************************************************************************/
618 {
619  // Return
620  return;
621 }
virtual std::string print(const GChatter &chatter=NORMAL) const
Print point source information.
virtual double eval(const double &offset, const bool &gradients=false) const
Evaluate function.
#define G_WRITE
virtual ~GCTAModelRadialGauss(void)
Destructor.
GCTAModelRadialGauss(void)
Void constructor.
const double & factor_gradient(void) const
Return parameter gradient factor.
double romberg(std::vector< double > bounds, const int &order=5)
Perform Romberg integration.
Definition: GIntegral.cpp:380
const std::string & name(void) const
Return parameter name.
virtual GCTAModelRadialGauss & operator=(const GCTAModelRadialGauss &model)
Assignment operator.
Random number generator class definition.
void free_members(void)
Delete class members.
const double pi
Definition: GMath.hpp:35
double gradient(void) const
Return parameter gradient.
void write(GXmlElement &xml) const
Set or update parameter attributes in XML element.
Definition: GModelPar.cpp:353
GVector cos(const GVector &vector)
Computes cosine of vector elements.
Definition: GVector.cpp:1100
XML element node class.
Definition: GXmlElement.hpp:47
void free_members(void)
Delete class members.
Random number generator class.
Definition: GRan.hpp:44
void init_members(void)
Initialise class members.
virtual int elements(void) const
Return number of GXMLElement children of node.
Definition: GXmlNode.cpp:580
virtual void clear(void)
Clear instance.
Gammalib tools definition.
virtual void read(const GXmlElement &xml)
Read model from XML element.
GIntegral class interface definition.
Definition: GIntegral.hpp:45
CTA radial model registry class definition.
virtual std::string type(void) const
Return model type.
double min(void) const
Return parameter minimum boundary.
Radial Gaussian CTA model class.
Interface definition for the spatial model registry class.
bool is_notanumber(const double &x)
Signal if argument is not a number.
Definition: GTools.hpp:184
const double & scale(void) const
Return parameter scale.
bool is_infinite(const double &x)
Signal if argument is infinite.
Definition: GTools.hpp:167
GModelPar m_sigma
Width parameter (degrees^2)
const double deg2rad
Definition: GMath.hpp:43
GVector sqrt(const GVector &vector)
Computes square root of vector elements.
Definition: GVector.cpp:1268
bool has_grad(void) const
Signal if parameter gradient is computed analytically.
#define G_READ
const GXmlAttribute * attribute(const int &index) const
Return attribute.
void free(void)
Free a parameter.
double uniform(void)
Returns random double precision floating value in range 0 to 1.
Definition: GRan.cpp:242
void clear(void)
Clear parameter.
CTA instrument direction class interface definition.
void copy_members(const GCTAModelRadialGauss &model)
Copy class members.
int size(void) const
Return number of model parameters.
GChatter
Definition: GTypemaps.hpp:33
double sigma(void) const
Return Gaussian width parameter.
virtual GCTAModelRadialGauss * clone(void) const
Clone instance.
CTA observation class interface definition.
virtual double mc_max_value(const GCTAObservation &obs) const
Return maximum function value for Monte Carlo simulations.
virtual GXmlElement * element(const int &index)
Return pointer to GXMLElement child.
Definition: GXmlNode.cpp:634
void read(const GXmlElement &xml)
Extract parameter attributes from XML element.
Definition: GModelPar.cpp:229
double value(void) const
Return parameter value.
bool has_min(void) const
Signal if parameter has minimum boundary.
const std::string & unit(void) const
Return parameter unit.
Interface definition for the CTA radial model registry class.
const GCTAModelRadialGauss g_cta_radial_gauss_seed
GVector sin(const GVector &vector)
Computes sine of vector elements.
Definition: GVector.cpp:1226
Exception handler interface definition.
Radial Gaussian model class interface definition.
void init_members(void)
Initialise class members.
GVector exp(const GVector &vector)
Computes exponential of vector elements.
Definition: GVector.cpp:1142
CTA instrument direction class.
Definition: GCTAInstDir.hpp:59
virtual GCTAInstDir mc(GRan &ran) const
Returns MC instrument direction.
Abstract radial acceptance model class.
virtual double omega(void) const
Returns integral over radial model (in steradians)
const double twopi
Definition: GMath.hpp:36
virtual GXmlNode * append(const GXmlNode &node)
Append XML child node.
Definition: GXmlNode.cpp:284
CTA observation class.
std::string parformat(const std::string &s, const int &indent=0)
Convert string in parameter format.
Definition: GTools.cpp:1022
std::vector< GModelPar * > m_pars
Parameter pointers.
Integration class interface definition.
virtual GCTAModelRadial & operator=(const GCTAModelRadial &model)
Assignment operator.
Spatial model registry class definition.
virtual void write(GXmlElement &xml) const
Write model into XML element.
Mathematical function definitions.
std::string str(const unsigned short int &value)
Convert unsigned short integer value into string.
Definition: GTools.cpp:413