GammaLib  2.0.0.dev
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
GCTAModelSpatialGradient.cpp
Go to the documentation of this file.
1 /***************************************************************************
2  * GCTAModelSpatialGradient.cpp - Spatial gradient CTA model class *
3  * ----------------------------------------------------------------------- *
4  * copyright (C) 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 GCTAModelSpatialGradient.cpp
23  * @brief Spatial gradient CTA interface definition
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 "GCTAObservation.hpp"
37 #include "GCTAInstDir.hpp"
40 
41 /* __ Constants __________________________________________________________ */
42 
43 /* __ Globals ____________________________________________________________ */
45 const GCTAModelSpatialRegistry g_cta_spatial_gradient_registry(&g_cta_spatial_gradient_seed);
46 
47 /* __ Method name definitions ____________________________________________ */
48 #define G_READ "GCTAModelSpatialGradient::read(GXmlElement&)"
49 #define G_WRITE "GCTAModelSpatialGradient::write(GXmlElement&)"
50 
51 /* __ Macros _____________________________________________________________ */
52 
53 /* __ Coding definitions _________________________________________________ */
54 //#define G_DEBUG_MC //!< Debug MC method
55 
56 /* __ Debug definitions __________________________________________________ */
57 
58 
59 /*==========================================================================
60  = =
61  = Constructors/destructors =
62  = =
63  ==========================================================================*/
64 
65 /***********************************************************************//**
66  * @brief Void constructor
67  ***************************************************************************/
69 {
70  // Initialise members
71  init_members();
72 
73  // Return
74  return;
75 }
76 
77 
78 /***********************************************************************//**
79  * @brief Gradient constructor
80  *
81  * @param[in] detx_gradient DETX gradient (degrees\f$^{-1}\f$).
82  * @param[in] dety_gradient DETY gradient (degrees\f$^{-1}\f$).
83  ***************************************************************************/
85  const double& dety_gradient) :
87 {
88  // Initialise members
89  init_members();
90 
91  // Assign gradients
92  this->detx_gradient(detx_gradient);
93  this->dety_gradient(dety_gradient);
94 
95  // Return
96  return;
97 }
98 
99 
100 /***********************************************************************//**
101  * @brief XML constructor
102  *
103  * @param[in] xml XML element.
104  *
105  * Creates instance of a spatial gradient model by extracting information
106  * from an XML element. See GCTAModelSpatialGradient::read() for more
107  * information about the expected structure of the XML element.
108  ***************************************************************************/
110 {
111  // Initialise members
112  init_members();
113 
114  // Read information from XML element
115  read(xml);
116 
117  // Return
118  return;
119 }
120 
121 
122 /***********************************************************************//**
123  * @brief Copy constructor
124  *
125  * @param[in] model Spatial gradient model.
126  ***************************************************************************/
128  GCTAModelSpatial(model)
129 {
130  // Initialise members
131  init_members();
132 
133  // Copy members
134  copy_members(model);
135 
136  // Return
137  return;
138 }
139 
140 
141 /***********************************************************************//**
142  * @brief Destructor
143  ***************************************************************************/
145 {
146  // Free members
147  free_members();
148 
149  // Return
150  return;
151 }
152 
153 
154 /*==========================================================================
155  = =
156  = Operators =
157  = =
158  ==========================================================================*/
159 
160 /***********************************************************************//**
161  * @brief Assignment operator
162  *
163  * @param[in] model Spatial gradient model.
164  ***************************************************************************/
166 {
167  // Execute only if object is not identical
168  if (this != &model) {
169 
170  // Copy base class members
171  this->GCTAModelSpatial::operator=(model);
172 
173  // Free members
174  free_members();
175 
176  // Initialise members
177  init_members();
178 
179  // Copy members
180  copy_members(model);
181 
182  } // endif: object was not identical
183 
184  // Return
185  return *this;
186 }
187 
188 
189 /*==========================================================================
190  = =
191  = Public methods =
192  = =
193  ==========================================================================*/
194 
195 /***********************************************************************//**
196  * @brief Clear instance
197  ***************************************************************************/
199 {
200  // Free class members (base and derived classes, derived class first)
201  free_members();
203 
204  // Initialise members
206  init_members();
207 
208  // Return
209  return;
210 }
211 
212 
213 /***********************************************************************//**
214  * @brief Clone instance
215  ***************************************************************************/
217 {
218  return new GCTAModelSpatialGradient(*this);
219 }
220 
221 
222 /***********************************************************************//**
223  * @brief Evaluate function
224  *
225  * @param[in] dir Event direction.
226  * @param[in] energy Event energy (not used).
227  * @param[in] time Event time (not used).
228  * @param[in] gradients Compute gradients?
229  * @return Function value
230  *
231  * Evaluates the spatial gradient model for a given event direction. The
232  * energy and time of the event are not used.
233  *
234  * The spatial gradient model is defined as
235  *
236  * \f[f(x,y) = 1 + g_x \times x + g_y \times y\f]
237  *
238  * where
239  * \f$x\f$ is x direction,
240  * \f$y\f$ is y direction,
241  * \f$g_x\f$ is the spatial gradient in the x direction, and
242  * \f$g_y\f$ is the spatial gradient in the y direction.
243  *
244  * If the @p gradients flag is true the method will also compute the partial
245  * derivatives of the parameters. The partial derivative of the spatial
246  * gradient model are given by
247  *
248  * \f[ \frac{df}{dg_{xv}} = g_{xs} x\f]
249  *
250  * and
251  *
252  * \f[ \frac{df}{dg_{yv}} = g_{ys} y\f]
253  *
254  * where
255  * \f$g_{xv}\f$ is the value part and \f$g_{xs}\f$ is the scaling part of
256  * gradient in x, and
257  * \f$g_{yv}\f$ is the value part and \f$g_{ys}\f$ is the scaling part of
258  * gradient in y.
259  ***************************************************************************/
261  const GEnergy& energy,
262  const GTime& time,
263  const bool& gradients) const
264 {
265  // Get detx and dety in degrees
266  double detx = dir.detx() * gammalib::rad2deg;
267  double dety = dir.dety() * gammalib::rad2deg;
268 
269  // Compute value
270  double value = 1.0 + m_detx_gradient.value() * detx +
271  m_dety_gradient.value() * dety;
272 
273  // Optionally compute partial derivatives
274  if (gradients) {
277  }
278 
279  // Compile option: Check for NaN/Inf
280  #if defined(G_NAN_CHECK)
281  if (gammalib::is_notanumber(value) || gammalib::is_infinite(value)) {
282  std::cout << "*** ERROR: GCTAModelSpatialGradient::eval";
283  std::cout << "(detx=" << detx << ", dety=" << dety;
284  std::cout << "): NaN/Inf encountered";
285  std::cout << " (value=" << value;
286  std::cout << ")" << std::endl;
287  }
288  #endif
289 
290  // Return value
291  return value;
292 }
293 
294 
295 /***********************************************************************//**
296  * @brief Return maximum function value for Monte Carlo simulations
297  *
298  * @param[in] obs CTA Observation.
299  * @return Maximum function value for Monte Carlo simulations.
300  ***************************************************************************/
302 {
303  // Get DETX and DETY value of RoI centre in degrees
304  GCTAInstDir roi_centre = obs.roi().centre();
305  double detx_centre = roi_centre.detx() * gammalib::rad2deg;
306  double dety_centre = roi_centre.dety() * gammalib::rad2deg;
307 
308  // Get DETX and DETY minima and maximum in degrees
309  double radius = obs.roi().radius();
310  double detx_min = detx_centre - radius;
311  double detx_max = detx_centre + radius;
312  double dety_min = dety_centre - radius;
313  double dety_max = dety_centre + radius;
314 
315  // Get maximum value
316  double value_min = m_detx_gradient.value() * detx_min;
317  double value_max = m_detx_gradient.value() * detx_max;
318  double valuex = (value_min > value_max) ? value_min : value_max;
319  value_min = m_dety_gradient.value() * dety_min;
320  value_max = m_dety_gradient.value() * dety_max;
321  double valuey = (value_min > value_max) ? value_min : value_max;
322  double value = 1.0 + valuex + valuey;
323 
324  // Return maximum value
325  return value;
326 }
327 
328 
329 /***********************************************************************//**
330  * @brief Read model from XML element
331  *
332  * @param[in] xml XML element.
333  *
334  * Read the gradient spatial model information from an XML element.
335  ***************************************************************************/
337 {
338  // Get parameter pointers
341 
342  // Read parameters
343  m_detx_gradient.read(*detx);
344  m_dety_gradient.read(*dety);
345 
346  // Return
347  return;
348 }
349 
350 
351 /***********************************************************************//**
352  * @brief Write model into XML element
353  *
354  * @param[in] xml XML element.
355  *
356  * @exception GException::invalid_value
357  * Spatial model is not of valid type.
358  *
359  * Write the gradient spatial model information into an XML element.
360  ***************************************************************************/
362 {
363  // Set model type
364  if (xml.attribute("type") == "") {
365  xml.attribute("type", type());
366  }
367 
368  // Verify model type
369  if (xml.attribute("type") != type()) {
370  std::string msg = "Spatial model \""+xml.attribute("type")+
371  "\" is not of type \""+type()+"\".";
373  }
374 
375  // Get XML parameters
378 
379  // Write parameters
380  m_detx_gradient.write(*detx);
381  m_dety_gradient.write(*dety);
382 
383  // Return
384  return;
385 }
386 
387 
388 /***********************************************************************//**
389  * @brief Print point source information
390  *
391  * @param[in] chatter Chattiness.
392  * @return String containing point source information.
393  ***************************************************************************/
394 std::string GCTAModelSpatialGradient::print(const GChatter& chatter) const
395 {
396  // Initialise result string
397  std::string result;
398 
399  // Continue only if chatter is not silent
400  if (chatter != SILENT) {
401 
402  // Append header
403  result.append("=== GCTAModelSpatialGradient ===");
404 
405  // Append information
406  result.append("\n"+gammalib::parformat("Number of parameters") +
407  gammalib::str(size()));
408  for (int i = 0; i < size(); ++i) {
409  result.append("\n"+m_pars[i]->print(chatter));
410  }
411 
412  } // endif: chatter was not silent
413 
414  // Return result
415  return result;
416 }
417 
418 
419 /*==========================================================================
420  = =
421  = Private methods =
422  = =
423  ==========================================================================*/
424 
425 /***********************************************************************//**
426  * @brief Initialise class members
427  ***************************************************************************/
429 {
430  // Initialise detx gradient
432  m_detx_gradient.name("Grad_DETX");
433  m_detx_gradient.unit("deg^-1");
434  m_detx_gradient.value(0.0);
436  m_detx_gradient.scale(1.0);
439 
440  // Initialise dety gradient
442  m_dety_gradient.name("Grad_DETY");
443  m_dety_gradient.unit("deg^-1");
444  m_dety_gradient.value(0.0);
446  m_dety_gradient.scale(1.0);
449 
450  // Set parameter pointer(s)
451  m_pars.clear();
452  m_pars.push_back(&m_detx_gradient);
453  m_pars.push_back(&m_dety_gradient);
454 
455  // Return
456  return;
457 }
458 
459 
460 /***********************************************************************//**
461  * @brief Copy class members
462  *
463  * @param[in] model Radial Gaussian model.
464  ***************************************************************************/
466 {
467  // Copy members
470 
471  // Set parameter pointer(s)
472  m_pars.clear();
473  m_pars.push_back(&m_detx_gradient);
474  m_pars.push_back(&m_dety_gradient);
475 
476  // Return
477  return;
478 }
479 
480 
481 /***********************************************************************//**
482  * @brief Delete class members
483  ***************************************************************************/
485 {
486  // Return
487  return;
488 }
GModelPar m_detx_gradient
DETX gradient.
Abstract spatial model class.
const double & factor_gradient(void) const
Return parameter factor gradient.
void detx(const double &x)
Set DETX coordinate (in radians)
virtual GCTAModelSpatialGradient * clone(void) const
Clone instance.
const std::string & name(void) const
Return parameter name.
Random number generator class definition.
void copy_members(const GCTAModelSpatialGradient &model)
Copy class members.
double gradient(void) const
Return parameter gradient.
#define G_READ
void write(GXmlElement &xml) const
Set or update parameter attributes in XML element.
Definition: GModelPar.cpp:351
virtual void clear(void)
Clear instance.
virtual GCTAModelSpatial & operator=(const GCTAModelSpatial &model)
Assignment operator.
const double & radius(void) const
Returns radius of region of interest in degrees.
Definition: GCTARoi.hpp:124
XML element node class.
Definition: GXmlElement.hpp:47
Time class.
Definition: GTime.hpp:54
Gammalib tools definition.
Interface definition for the spatial model registry class.
const GCTAModelSpatialGradient g_cta_spatial_gradient_seed
void free_members(void)
Delete class members.
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
bool has_grad(void) const
Signal if parameter gradient is computed analytically.
virtual void read(const GXmlElement &xml)
Read model from XML element.
const GXmlAttribute * attribute(const int &index) const
Return attribute.
GCTARoi roi(void) const
Get Region of Interest.
#define G_WRITE
void free(void)
Free a parameter.
GModelPar m_dety_gradient
DETY gradient.
const GCTAInstDir & centre(void) const
Returns region of interest centre.
Definition: GCTARoi.hpp:110
void free_members(void)
Delete class members.
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:1527
GCTAModelSpatialGradient(void)
Void constructor.
virtual ~GCTAModelSpatialGradient(void)
Destructor.
void init_members(void)
Initialise class members.
void clear(void)
Clear parameter.
CTA instrument direction class interface definition.
int size(void) const
Return number of model parameters.
GChatter
Definition: GTypemaps.hpp:33
Spatial gradient CTA interface definition.
CTA observation class interface definition.
virtual void write(GXmlElement &xml) const
Write model into XML element.
virtual double eval(const GCTAInstDir &dir, const GEnergy &energy, const GTime &time, const bool &gradients=false) const
Evaluate function.
virtual double mc_max_value(const GCTAObservation &obs) const
Return maximum function value for Monte Carlo simulations.
void init_members(void)
Initialise class members.
Spatial gradient CTA model class.
double detx_gradient(void) const
Return DETX gradient.
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 dety_gradient(void) const
Return DETY gradient.
virtual GCTAModelSpatialGradient & operator=(const GCTAModelSpatialGradient &model)
Assignment operator.
CTA instrument direction class.
Definition: GCTAInstDir.hpp:63
void dety(const double &y)
Set DETY coordinate (in radians)
const double rad2deg
Definition: GMath.hpp:44
CTA observation class.
std::string parformat(const std::string &s, const int &indent=0)
Convert string in parameter format.
Definition: GTools.cpp:1033
std::vector< GModelPar * > m_pars
Parameter pointers.
virtual std::string type(void) const
Return model type.
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:1579
Spatial model registry class definition.
Mathematical function definitions.
virtual std::string print(const GChatter &chatter=NORMAL) const
Print point source information.
Class that handles energies in a unit independent way.
Definition: GEnergy.hpp:48
std::string str(const unsigned short int &value)
Convert unsigned short integer value into string.
Definition: GTools.cpp:415