GammaLib  1.7.0.dev
 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-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 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  * @exception GException::model_invalid_parnum
413  * Invalid number of model parameters found in XML element.
414  * @exception GException::model_invalid_parvalue
415  * Non-positive parameter value found.
416  * @exception GException::model_invalid_parlimit
417  * Missing or non-positive minimum parameter boundary.
418  * @exception GException::model_invalid_parnames
419  * Invalid model parameter names found in XML element.
420  *
421  * Read the Profile radial model information from an XML element. The XML
422  * element is required to have 3 parameter named "Width", "Core", and "Tail".
423  ***************************************************************************/
425 {
426  // Verify that XML element has exactly 3 parameter
427  if (xml.elements() != 3 || xml.elements("parameter") != 3) {
429  "Radial Profile model requires exactly 3 parameters.");
430  }
431 
432  // Extract model parameters
433  int npar[] = {0, 0, 0};
434  for (int i = 0; i < 3; ++i) {
435 
436  // Get parameter element
437  const GXmlElement* par = xml.element("parameter", i);
438 
439  // Handle Width
440  if (par->attribute("name") == "Width") {
441  m_width.read(*par);
442  if (m_width.value() <= 0.0) {
444  "\"Width\" parameter is required to be positive.");
445  }
446  if (!m_width.has_min() || m_width.min() <= 0.0) {
448  "\"Width\" parameter requires positive minimum boundary.");
449  }
450  npar[0]++;
451  }
452 
453  // Handle Core
454  if (par->attribute("name") == "Core") {
455  m_core.read(*par);
456  if (m_core.value() <= 0.0) {
458  "\"Core\" parameter is required to be positive.");
459  }
460  if (!m_core.has_min() || m_core.min() <= 0.0) {
462  "\"Core\" parameter requires positive minimum boundary.");
463  }
464  npar[1]++;
465  }
466 
467  // Handle Tail
468  if (par->attribute("name") == "Tail") {
469  m_tail.read(*par);
470  if (m_tail.value() <= 0.0) {
472  "\"Tail\" parameter is required to be positive.");
473  }
474  if (!m_tail.has_min() || m_tail.min() <= 0.0) {
476  "\"Tail\" parameter requires positive minimum boundary.");
477  }
478  npar[2]++;
479  }
480 
481  } // endfor: looped over all parameters
482 
483  // Verify that all parameters were found
484  if (npar[0] != 1 || npar[1] != 1 || npar[2] != 1) {
486  "Require \"Width\", \"Core\" and \"Tail\" parameters.");
487  }
488 
489  // Return
490  return;
491 }
492 
493 
494 /***********************************************************************//**
495  * @brief Write model into XML element
496  *
497  * @param[in] xml XML element.
498  *
499  * @exception GException::model_invalid_spatial
500  * Existing XML element is not of type 'ProfileFunction'
501  * @exception GException::model_invalid_parnum
502  * Invalid number of model parameters found in XML element.
503  * @exception GException::model_invalid_parnames
504  * Invalid model parameter names found in XML element.
505  *
506  * Write the radial Profile model information into an XML element. The XML
507  * element will have 3 parameter leafs named "Width", "Core" and "Tail".
508  ***************************************************************************/
510 {
511  // Set model type
512  if (xml.attribute("type") == "") {
513  xml.attribute("type", type());
514  }
515 
516  // Verify model type
517  if (xml.attribute("type") != type()) {
519  "Radial Profile model is not of type \""+type()+"\".");
520  }
521 
522  // If XML element has 0 nodes then append 3 parameter nodes
523  if (xml.elements() == 0) {
524  xml.append(GXmlElement("parameter name=\"Width\""));
525  xml.append(GXmlElement("parameter name=\"Core\""));
526  xml.append(GXmlElement("parameter name=\"Tail\""));
527  }
528 
529  // Verify that XML element has exactly 3 parameters
530  if (xml.elements() != 3 || xml.elements("parameter") != 3) {
532  "Radial Profile model requires exactly 3 parameters.");
533  }
534 
535  // Set or update model parameter attributes
536  int npar[] = {0, 0, 0};
537  for (int i = 0; i < 3; ++i) {
538 
539  // Get parameter element
540  GXmlElement* par = xml.element("parameter", i);
541 
542  // Handle prefactor
543  if (par->attribute("name") == "Width") {
544  m_width.write(*par);
545  npar[0]++;
546  }
547 
548  // Handle index
549  else if (par->attribute("name") == "Core") {
550  m_core.write(*par);
551  npar[1]++;
552  }
553 
554  // Handle pivot energy
555  else if (par->attribute("name") == "Tail") {
556  m_tail.write(*par);
557  npar[2]++;
558  }
559 
560  } // endfor: looped over all parameters
561 
562  // Check of all required parameters are present
563  if (npar[0] != 1 || npar[1] != 1 || npar[2] != 1) {
565  "Require \"Width\", \"Core\" and \"Tail\" parameters.");
566  }
567 
568  // Return
569  return;
570 }
571 
572 
573 /***********************************************************************//**
574  * @brief Print point source information
575  *
576  * @param[in] chatter Chattiness (defaults to NORMAL).
577  * @return String containing point source information.
578  ***************************************************************************/
579 std::string GCTAModelRadialProfile::print(const GChatter& chatter) const
580 {
581  // Initialise result string
582  std::string result;
583 
584  // Continue only if chatter is not silent
585  if (chatter != SILENT) {
586 
587  // Append header
588  result.append("=== GCTAModelRadialProfile ===");
589 
590  // Append information
591  result.append("\n"+gammalib::parformat("Number of parameters") +
592  gammalib::str(size()));
593  for (int i = 0; i < size(); ++i) {
594  result.append("\n"+m_pars[i]->print(chatter));
595  }
596 
597  } // endif: chatter was not silent
598 
599  // Return result
600  return result;
601 }
602 
603 
604 /*==========================================================================
605  = =
606  = Private methods =
607  = =
608  ==========================================================================*/
609 
610 /***********************************************************************//**
611  * @brief Initialise class members
612  ***************************************************************************/
614 {
615  // Initialise width parameter
616  m_width.clear();
617  m_width.name("Width");
618  m_width.unit("deg");
619  m_width.scale(1.0);
620  m_width.value(1.5); // default: 1.5 deg
621  m_width.min(0.1); // min: 0.1 deg
622  m_width.free();
623  m_width.gradient(0.0);
624  m_width.has_grad(false);
625 
626  // Initialise core parameter
627  m_core.clear();
628  m_core.name("Core");
629  m_core.scale(1.0);
630  m_core.value(3.0); // default: 3.0
631  m_core.min(1.0); // min: 1.0 (could even be larger)
632  m_core.free();
633  m_core.gradient(0.0);
634  m_core.has_grad(false);
635 
636  // Initialise tail parameter
637  m_tail.clear();
638  m_tail.name("Tail");
639  m_tail.scale(1.0);
640  m_tail.value(5.0); // default: 5.0
641  m_tail.min(1.0); // min: 1.0
642  m_tail.fix();
643  m_tail.gradient(0.0);
644  m_tail.has_grad(false);
645 
646  // Set parameter pointer(s)
647  m_pars.clear();
648  m_pars.push_back(&m_width);
649  m_pars.push_back(&m_core);
650  m_pars.push_back(&m_tail);
651 
652  // Return
653  return;
654 }
655 
656 
657 /***********************************************************************//**
658  * @brief Copy class members
659  *
660  * @param[in] model Radial Profileian model.
661  ***************************************************************************/
663 {
664  // Copy members
665  m_width = model.m_width;
666  m_core = model.m_core;
667  m_tail = model.m_tail;
668 
669  // Set parameter pointer(s)
670  m_pars.clear();
671  m_pars.push_back(&m_width);
672  m_pars.push_back(&m_core);
673  m_pars.push_back(&m_tail);
674 
675  // Return
676  return;
677 }
678 
679 
680 /***********************************************************************//**
681  * @brief Delete class members
682  ***************************************************************************/
684 {
685  // Return
686  return;
687 }
688 
689 
690 /*==========================================================================
691  = =
692  = Friends =
693  = =
694  ==========================================================================*/
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:380
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:353
GVector cos(const GVector &vector)
Computes cosine of vector elements.
Definition: GVector.cpp:1100
double tail(void) const
Return profile tail.
#define G_READ
XML element node class.
Definition: GXmlElement.hpp:47
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.
virtual int elements(void) const
Return number of GXMLElement children of node.
Definition: GXmlNode.cpp:580
Gammalib tools definition.
GIntegral class interface definition.
Definition: GIntegral.hpp:45
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:184
const double & scale(void) const
Return parameter scale.
bool is_infinite(const double &x)
Signal if argument is infinite.
Definition: GTools.hpp:167
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.
const GXmlAttribute * attribute(const int &index) const
Return attribute.
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.
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.
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 GXmlElement * element(const int &index)
Return pointer to GXMLElement child.
Definition: GXmlNode.cpp:634
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.
bool has_min(void) const
Signal if parameter has minimum boundary.
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:1332
Interface definition for the CTA radial model registry class.
GVector sin(const GVector &vector)
Computes sine of vector elements.
Definition: GVector.cpp:1226
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:59
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.
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.
virtual std::string type(void) const
Return model type.
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:413