GammaLib  2.0.0
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
GCTAModelRadialProfile.cpp
Go to the documentation of this file.
1 /***************************************************************************
2  * GCTAModelRadialProfile.cpp - Radial Profile CTA model class *
3  * ----------------------------------------------------------------------- *
4  * copyright (C) 2011-2021 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 GCTAModelRadialProfile.cpp
23  * @brief Radial Profile 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"
42 
43 /* __ Constants __________________________________________________________ */
44 
45 /* __ Globals ____________________________________________________________ */
47 const GCTAModelRadialRegistry g_cta_radial_profile_registry(&g_cta_radial_profile_seed);
48 const GCTAModelSpatialRegistry g_cta_radial_profile_spatial_registry(&g_cta_radial_profile_seed);
49 
50 /* __ Method name definitions ____________________________________________ */
51 #define G_READ "GCTAModelRadialProfile::read(GXmlElement&)"
52 #define G_WRITE "GCTAModelRadialProfile::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] width Width.
85  * @param[in] core Core size.
86  * @param[in] tail Tail size.
87  ***************************************************************************/
89  const double& core,
90  const double& tail) : GCTAModelRadial()
91 {
92  // Initialise members
93  init_members();
94 
95  // Assign values
96  this->width(width);
97  this->core(core);
98  this->tail(tail);
99 
100  // Return
101  return;
102 }
103 
104 
105 /***********************************************************************//**
106  * @brief XML constructor
107  *
108  * @param[in] xml XML element.
109  *
110  * Creates instance of a radial Profile model by extracting information from
111  * an XML element. See GCTAModelRadialProfile::read() for more information
112  * about the expected structure of the XML element.
113  ***************************************************************************/
115  : GCTAModelRadial()
116 {
117  // Initialise members
118  init_members();
119 
120  // Read information from XML element
121  read(xml);
122 
123  // Return
124  return;
125 }
126 
127 
128 /***********************************************************************//**
129  * @brief Copy constructor
130  *
131  * @param[in] model Radial Profile model.
132  ***************************************************************************/
134  : GCTAModelRadial(model)
135 {
136  // Initialise members
137  init_members();
138 
139  // Copy members
140  copy_members(model);
141 
142  // Return
143  return;
144 }
145 
146 
147 /***********************************************************************//**
148  * @brief Destructor
149  ***************************************************************************/
151 {
152  // Free members
153  free_members();
154 
155  // Return
156  return;
157 }
158 
159 
160 /*==========================================================================
161  = =
162  = Operators =
163  = =
164  ==========================================================================*/
165 
166 /***********************************************************************//**
167  * @brief Assignment operator
168  *
169  * @param[in] model Radial Profile model.
170  ***************************************************************************/
172 {
173  // Execute only if object is not identical
174  if (this != &model) {
175 
176  // Copy base class members
177  this->GCTAModelRadial::operator=(model);
178 
179  // Free members
180  free_members();
181 
182  // Initialise members
183  init_members();
184 
185  // Copy members
186  copy_members(model);
187 
188  } // endif: object was not identical
189 
190  // Return
191  return *this;
192 }
193 
194 
195 /*==========================================================================
196  = =
197  = Public methods =
198  = =
199  ==========================================================================*/
200 
201 /***********************************************************************//**
202  * @brief Clear instance
203 ***************************************************************************/
205 {
206  // Free class members (base and derived classes, derived class first)
207  free_members();
209 
210  // Initialise members
212  init_members();
213 
214  // Return
215  return;
216 }
217 
218 
219 /***********************************************************************//**
220  * @brief Clone instance
221 ***************************************************************************/
223 {
224  return new GCTAModelRadialProfile(*this);
225 }
226 
227 
228 /***********************************************************************//**
229  * @brief Evaluate function
230  *
231  * @param[in] offset Offset angle [degrees].
232  * @param[in] gradients Compute gradients?
233  * @return Function value
234  *
235  * Evaluates the Profile model for a given offset. The Profile model is
236  * defined as
237  * \f[f(\theta) = (1 + (\theta/c_0)^{c_1})^{-c_2/c_1}\f]
238  * where
239  * \f$\theta\f$ is the offset angle (in degrees),
240  * \f$c_0\f$ is the width of the profile (width),
241  * \f$c_1\f$ is the width of the central plateau (core), and
242  * \f$c_2\f$ is the size of the tail (tail).
243  *
244  * Note that no analytical partial derivatives are implemented for this
245  * function.
246  ***************************************************************************/
247 double GCTAModelRadialProfile::eval(const double& offset,
248  const bool& gradients) const
249 {
250  // Compute value
251  double arg = 1.0 + std::pow(offset / width(), core());
252  double value = std::pow(arg, -tail()/core());
253 
254  // Compile option: Check for NaN/Inf
255  #if defined(G_NAN_CHECK)
256  if (gammalib::is_notanumber(value) || gammalib::is_infinite(value)) {
257  std::cout << "*** ERROR: GCTAModelRadialProfile::eval";
258  std::cout << "(offset=" << offset << "): NaN/Inf encountered";
259  std::cout << " (value=" << value;
260  std::cout << ", width=" << width();
261  std::cout << ", core=" << core();
262  std::cout << ", tail=" << tail();
263  std::cout << ")" << std::endl;
264  }
265  #endif
266 
267  // Return value
268  return value;
269 }
270 
271 
272 /***********************************************************************//**
273  * @brief Returns MC instrument direction
274  *
275  * @param[in,out] ran Random number generator.
276  * @return CTA instrument direction.
277  *
278  * Draws an arbitrary CTA instrument position from
279  * \f[f(\theta) = \sin(\theta) (1 + (\theta/c_0)^{c_1})^{-c_2/c_1}\f]
280  * where
281  * \f$\theta\f$ is the offset angle (in degrees),
282  * \f$c_0\f$ is the width of the profile (width),
283  * \f$c_1\f$ is the width of the central plateau (core), and
284  * \f$c_2\f$ is the size of the tail (tail).
285  * using the rejection method.
286  *
287  * The maximum offset angle that is thrown by this method is fixed to 10
288  * degrees.
289  *
290  * @todo Method can be optimised by using a random deviate of sin instead
291  * of the uniform random deviate which leads to many unnecessary
292  * rejections.
293  ***************************************************************************/
295 {
296  // Set constants
297  const double offset_max = 10.0;
298  const double u_max = std::sin(offset_max * gammalib::deg2rad);
299 
300  // Debug option: initialise number if samples thrown for one value
301  #if defined(G_DEBUG_MC)
302  int n_samples = 0;
303  #endif
304 
305  // Simulate offset from photon arrival direction until we're not
306  // rejected anymore
307  double value = 0.0;
308  double u = 1.0;
309  double offset = 0.0;
310  do {
311  // Throw offset angle
312  offset = ran.uniform() * offset_max;
313 
314  // Compute function value at this offset angle
315  value = std::sin(offset * gammalib::deg2rad) * eval(offset);
316 
317  // Throw value for rejection method
318  u = ran.uniform() * u_max;
319 
320  // Debug option: update number of samples
321  #if defined(G_DEBUG_MC)
322  n_samples++;
323  #endif
324 
325  } while (u > value);
326 
327  // Debug option: print number if samples thrown for one value
328  #if defined(G_DEBUG_MC)
329  std::cout << "#=" << n_samples << " ";
330  #endif
331 
332  // Simulate azimuth angle
333  double phi = 360.0 * ran.uniform();
334 
335  // Convert from degrees to radians
336  offset *= gammalib::deg2rad;
337  phi *= gammalib::deg2rad;
338 
339  // Compute DETX and DETY coordinates
340  double detx(0.0);
341  double dety(0.0);
342  if (offset > 0.0 ) {
343  detx = offset * std::cos(phi);
344  dety = offset * std::sin(phi);
345  }
346 
347  // Set instrument direction
348  GCTAInstDir dir(detx, dety);
349 
350  // Return instrument direction
351  return dir;
352 }
353 
354 
355 /***********************************************************************//**
356  * @brief Return maximum function value for Monte Carlo simulations
357  *
358  * @param[in] obs CTA Observation (not used).
359  * @return Maximum function value for Monte Carlo simulations.
360  ***************************************************************************/
362 {
363  // Set constants
364  const double offset_max = 10.0;
365  const double u_max = std::sin(offset_max * gammalib::deg2rad);
366 
367  // Return maximum value
368  return u_max;
369 }
370 
371 
372 /***********************************************************************//**
373  * @brief Returns integral over radial model (in steradians)
374  *
375  * Computes
376  * \f[\Omega = 2 \pi \int_0^{\pi} \sin \theta f(\theta) d\theta\f]
377  * where
378  * \f[f(\theta) = (1 + (\theta/c_0)^{c_1})^{-c_2/c_1}\f]
379  * \f$\theta\f$ is the offset angle (in degrees),
380  * \f$c_0\f$ is the width of the profile (width),
381  * \f$c_1\f$ is the width of the central plateau (core), and
382  * \f$c_2\f$ is the size of the tail (tail).
383  *
384  * The integration is performed numerically, and the upper integration bound
385  * \f$\pi\f$
386  * is fixed to 10 degrees.
387  ***************************************************************************/
389 {
390  // Set constants
391  const double offset_max_rad = 10.0 * gammalib::deg2rad;
392 
393  // Allocate integrand
395 
396  // Allocate intergal
397  GIntegral integral(&integrand);
398 
399  // Perform numerical integration
400  double omega = integral.romberg(0.0, offset_max_rad) * gammalib::twopi;
401 
402  // Return integral
403  return omega;
404 }
405 
406 
407 /***********************************************************************//**
408  * @brief Read model from XML element
409  *
410  * @param[in] xml XML element.
411  *
412  * Read the Profile radial model information from an XML element in the
413  * following format
414  *
415  * <radialModel type="...">
416  * <parameter name="Width" scale="1.0" value="1.5" min="0.1" max="10.0" free="1"/>
417  * <parameter name="Core" scale="1.0" value="3.0" min="1.0" max="10.0" free="1"/>
418  * <parameter name="Tail" scale="1.0" value="5.0" min="1.0" max="10.0" free="1"/>
419  * </radialModel>
420  ***************************************************************************/
422 {
423  // Verify that XML element has exactly 3 parameters
425 
426  // Get parameters
430 
431  // Read parameters
432  m_width.read(*width);
433  m_core.read(*core);
434  m_tail.read(*tail);
435 
436  // Return
437  return;
438 }
439 
440 
441 /***********************************************************************//**
442  * @brief Write model into XML element
443  *
444  * @param[in] xml XML element.
445  *
446  * Write the radial Profile model information into an XML element in the
447  * following format
448  *
449  * <radialModel type="...">
450  * <parameter name="Width" scale="1.0" value="1.5" min="0.1" max="10.0" free="1"/>
451  * <parameter name="Core" scale="1.0" value="3.0" min="1.0" max="10.0" free="1"/>
452  * <parameter name="Tail" scale="1.0" value="5.0" min="1.0" max="10.0" free="1"/>
453  * </radialModel>
454  ***************************************************************************/
456 {
457  // Check model type
459 
460  // Get or create parameters
464 
465  // Write parameters
466  m_width.write(*width);
467  m_core.write(*core);
468  m_tail.write(*tail);
469 
470  // Return
471  return;
472 }
473 
474 
475 /***********************************************************************//**
476  * @brief Print point source information
477  *
478  * @param[in] chatter Chattiness.
479  * @return String containing point source information.
480  ***************************************************************************/
481 std::string GCTAModelRadialProfile::print(const GChatter& chatter) const
482 {
483  // Initialise result string
484  std::string result;
485 
486  // Continue only if chatter is not silent
487  if (chatter != SILENT) {
488 
489  // Append header
490  result.append("=== GCTAModelRadialProfile ===");
491 
492  // Append information
493  result.append("\n"+gammalib::parformat("Number of parameters") +
494  gammalib::str(size()));
495  for (int i = 0; i < size(); ++i) {
496  result.append("\n"+m_pars[i]->print(chatter));
497  }
498 
499  } // endif: chatter was not silent
500 
501  // Return result
502  return result;
503 }
504 
505 
506 /*==========================================================================
507  = =
508  = Private methods =
509  = =
510  ==========================================================================*/
511 
512 /***********************************************************************//**
513  * @brief Initialise class members
514  ***************************************************************************/
516 {
517  // Initialise width parameter
518  m_width.clear();
519  m_width.name("Width");
520  m_width.unit("deg");
521  m_width.scale(1.0);
522  m_width.value(1.5); // default: 1.5 deg
523  m_width.min(0.1); // min: 0.1 deg
524  m_width.free();
525  m_width.gradient(0.0);
526  m_width.has_grad(false);
527 
528  // Initialise core parameter
529  m_core.clear();
530  m_core.name("Core");
531  m_core.scale(1.0);
532  m_core.value(3.0); // default: 3.0
533  m_core.min(1.0); // min: 1.0 (could even be larger)
534  m_core.free();
535  m_core.gradient(0.0);
536  m_core.has_grad(false);
537 
538  // Initialise tail parameter
539  m_tail.clear();
540  m_tail.name("Tail");
541  m_tail.scale(1.0);
542  m_tail.value(5.0); // default: 5.0
543  m_tail.min(1.0); // min: 1.0
544  m_tail.fix();
545  m_tail.gradient(0.0);
546  m_tail.has_grad(false);
547 
548  // Set parameter pointer(s)
549  m_pars.clear();
550  m_pars.push_back(&m_width);
551  m_pars.push_back(&m_core);
552  m_pars.push_back(&m_tail);
553 
554  // Return
555  return;
556 }
557 
558 
559 /***********************************************************************//**
560  * @brief Copy class members
561  *
562  * @param[in] model Radial Profileian model.
563  ***************************************************************************/
565 {
566  // Copy members
567  m_width = model.m_width;
568  m_core = model.m_core;
569  m_tail = model.m_tail;
570 
571  // Set parameter pointer(s)
572  m_pars.clear();
573  m_pars.push_back(&m_width);
574  m_pars.push_back(&m_core);
575  m_pars.push_back(&m_tail);
576 
577  // Return
578  return;
579 }
580 
581 
582 /***********************************************************************//**
583  * @brief Delete class members
584  ***************************************************************************/
586 {
587  // Return
588  return;
589 }
590 
591 
592 /*==========================================================================
593  = =
594  = Friends =
595  = =
596  ==========================================================================*/
virtual std::string print(const GChatter &chatter=NORMAL) const
Print point source information.
double romberg(std::vector< double > bounds, const int &order=5)
Perform Romberg integration.
Definition: GIntegral.cpp:381
const std::string & name(void) const
Return parameter name.
Random number generator class definition.
GCTAModelRadialProfile(void)
Void constructor.
void free_members(void)
Delete class members.
double gradient(void) const
Return parameter gradient.
void init_members(void)
Initialise class members.
void write(GXmlElement &xml) const
Set or update parameter attributes in XML element.
Definition: GModelPar.cpp:351
GVector cos(const GVector &vector)
Computes cosine of vector elements.
Definition: GVector.cpp:1190
double tail(void) const
Return profile tail.
#define G_READ
XML element node class.
Definition: GXmlElement.hpp:48
virtual void clear(void)
Clear instance.
virtual GCTAModelRadialProfile * clone(void) const
Clone instance.
Random number generator class.
Definition: GRan.hpp:44
void init_members(void)
Initialise class members.
Gammalib tools definition.
GIntegral class interface definition.
Definition: GIntegral.hpp:46
CTA radial model registry class definition.
double width(void) const
Return profile width.
double min(void) const
Return parameter minimum boundary.
Interface definition for the spatial model registry class.
bool is_notanumber(const double &x)
Signal if argument is not a number.
Definition: GTools.hpp:201
const double & scale(void) const
Return parameter scale.
bool is_infinite(const double &x)
Signal if argument is infinite.
Definition: GTools.hpp:184
virtual double omega(void) const
Returns integral over radial model (in steradians)
const double deg2rad
Definition: GMath.hpp:43
bool has_grad(void) const
Signal if parameter gradient is computed analytically.
void free_members(void)
Delete class members.
void free(void)
Free a parameter.
virtual double eval(const double &offset, const bool &gradients=false) const
Evaluate function.
void fix(void)
Fix a parameter.
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:1637
virtual ~GCTAModelRadialProfile(void)
Destructor.
double uniform(void)
Returns random double precision floating value in range 0 to 1.
Definition: GRan.cpp:242
void clear(void)
Clear parameter.
GModelPar m_core
Core parameter.
CTA instrument direction class interface definition.
int size(void) const
Return number of model parameters.
GChatter
Definition: GTypemaps.hpp:33
CTA observation class interface definition.
void xml_check_parnum(const std::string &origin, const GXmlElement &xml, const int &number)
Checks number of parameters.
Definition: GTools.cpp:1777
Radial Profile CTA model class.
virtual double mc_max_value(const GCTAObservation &obs) const
Return maximum function value for Monte Carlo simulations.
GModelPar m_tail
Tail parameter.
virtual GCTAInstDir mc(GRan &ran) const
Returns MC instrument direction.
virtual GCTAModelRadialProfile & operator=(const GCTAModelRadialProfile &model)
Assignment operator.
#define G_WRITE
const GCTAModelRadialProfile g_cta_radial_profile_seed
virtual void read(const GXmlElement &xml)
Read model from XML element.
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.
GVector pow(const GVector &vector, const double &power)
Computes tanh of vector elements.
Definition: GVector.cpp:1422
Interface definition for the CTA radial model registry class.
GVector sin(const GVector &vector)
Computes sine of vector elements.
Definition: GVector.cpp:1316
virtual void write(GXmlElement &xml) const
Write model into XML element.
Exception handler interface definition.
Radial Profile model class interface definition.
CTA instrument direction class.
Definition: GCTAInstDir.hpp:63
GModelPar m_width
Width parameter.
Abstract radial acceptance model class.
const double twopi
Definition: GMath.hpp:36
void copy_members(const GCTAModelRadialProfile &model)
Copy class members.
double core(void) const
Return profile core.
CTA observation class.
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.
Integration class interface definition.
virtual GCTAModelRadial & operator=(const GCTAModelRadial &model)
Assignment operator.
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:1689
Spatial model registry class definition.
Mathematical function definitions.
std::string str(const unsigned short int &value)
Convert unsigned short integer value into string.
Definition: GTools.cpp:489
void xml_check_type(const std::string &origin, GXmlElement &xml, const std::string &type)
Checks the model type.
Definition: GTools.cpp:1819