GammaLib  1.7.0.dev
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
GCTAModelSpatial.cpp
Go to the documentation of this file.
1 /***************************************************************************
2  * GCTAModelSpatial.cpp - Spatial model abstract base class *
3  * ----------------------------------------------------------------------- *
4  * copyright (C) 2018-2020 by Jurgen Knodlseder *
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 GCTAModelSpatial.cpp
23  * @brief Abstract spatial model class implementation
24  * @author Juergen Knoedlseder
25  */
26 
27 /* __ Includes ___________________________________________________________ */
28 #ifdef HAVE_CONFIG_H
29 #include <config.h>
30 #endif
31 #include "GException.hpp"
32 #include "GTools.hpp"
33 #include "GIntegral.hpp"
34 #include "GObservation.hpp"
35 #include "GCTASupport.hpp"
36 #include "GCTAEventList.hpp"
37 #include "GCTAInstDir.hpp"
38 #include "GCTAModelSpatial.hpp"
39 #include "GCTAResponse_helpers.hpp"
40 
41 /* __ Method name definitions ____________________________________________ */
42 #define G_ACCESS1 "GCTAModelSpatial::operator[](int&)"
43 #define G_ACCESS2 "GCTAModelSpatial::operator[](std::string&)"
44 #define G_MC "GCTAModelSpatial::mc(GEnergy&, GTime&, GCTAObservation&, "\
45  "GRan&)"
46 #define G_NPRED "GCTAModelSpatial::npred(GEnergy&, GTime&, GObservation&)"
47 
48 /* __ Macros _____________________________________________________________ */
49 
50 /* __ Coding definitions _________________________________________________ */
51 
52 /* __ Debug definitions __________________________________________________ */
53 
54 /* __ Constants __________________________________________________________ */
55 const double g_npred_resolution = 0.1*gammalib::deg2rad; //!< Scale of bkg.
56  // variation
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 Copy constructor
80  *
81  * @param[in] model Spatial model.
82  ***************************************************************************/
84 {
85  // Initialise members
86  init_members();
87 
88  // Copy members
89  copy_members(model);
90 
91  // Return
92  return;
93 }
94 
95 
96 /***********************************************************************//**
97  * @brief Destructor
98  ***************************************************************************/
100 {
101  // Free members
102  free_members();
103 
104  // Return
105  return;
106 }
107 
108 
109 /*==========================================================================
110  = =
111  = Operators =
112  = =
113  ==========================================================================*/
114 
115 /***********************************************************************//**
116  * @brief Assignment operator
117  *
118  * @param[in] model Spatial model.
119  ***************************************************************************/
121 {
122  // Execute only if object is not identical
123  if (this != &model) {
124 
125  // Free members
126  free_members();
127 
128  // Initialise members
129  init_members();
130 
131  // Copy members
132  copy_members(model);
133 
134  } // endif: object was not identical
135 
136  // Return
137  return *this;
138 }
139 
140 
141 /***********************************************************************//**
142  * @brief Returns model parameter
143  *
144  * @param[in] index Parameter index [0,...,size()-1].
145  *
146  * @exception GException::out_of_range
147  * Parameter index is out of range.
148  ***************************************************************************/
150 {
151  // Compile option: raise exception if index is out of range
152  #if defined(G_RANGE_CHECK)
153  if (index < 0 || index >= size()) {
154  throw GException::out_of_range(G_ACCESS1, "Spatial parameter index",
155  index, size());
156  }
157  #endif
158 
159  // Return reference
160  return *(m_pars[index]);
161 }
162 
163 
164 /***********************************************************************//**
165  * @brief Returns model parameter (const version)
166  *
167  * @param[in] index Parameter index [0,...,size()-1].
168  *
169  * @exception GException::out_of_range
170  * Parameter index is out of range.
171  ***************************************************************************/
172 const GModelPar& GCTAModelSpatial::operator[](const int& index) const
173 {
174  // Compile option: raise exception if index is out of range
175  #if defined(G_RANGE_CHECK)
176  if (index < 0 || index >= size()) {
177  throw GException::out_of_range(G_ACCESS1, "Spatial parameter index",
178  index, size());
179  }
180  #endif
181 
182  // Return reference
183  return *(m_pars[index]);
184 }
185 
186 
187 /***********************************************************************//**
188  * @brief Returns reference to model parameter
189  *
190  * @param[in] name Parameter name.
191  *
192  * @exception GException::invalid_argument
193  * Parameter with specified name not found in container.
194  ***************************************************************************/
195 GModelPar& GCTAModelSpatial::operator[](const std::string& name)
196 {
197  // Get parameter index
198  int index = 0;
199  for (; index < size(); ++index) {
200  if (m_pars[index]->name() == name)
201  break;
202  }
203 
204  // Throw exception if parameter name was not found
205  if (index >= size()) {
206  std::string msg = "Parameter \""+name+"\" not found in spatial "
207  "component of background model.";
209  }
210 
211  // Return reference
212  return *(m_pars[index]);
213 }
214 
215 
216 /***********************************************************************//**
217  * @brief Returns reference to model parameter (const version)
218  *
219  * @param[in] name Parameter name.
220  *
221  * @exception GException::invalid_argument
222  * Parameter with specified name not found in container.
223  ***************************************************************************/
224 const GModelPar& GCTAModelSpatial::operator[](const std::string& name) const
225 {
226  // Get parameter index
227  int index = 0;
228  for (; index < size(); ++index) {
229  if (m_pars[index]->name() == name)
230  break;
231  }
232 
233  // Throw exception if parameter name was not found
234  if (index >= size()) {
235  std::string msg = "Parameter \""+name+"\" not found in spatial "
236  "component of background model.";
238  }
239 
240  // Return reference
241  return *(m_pars[index]);
242 }
243 
244 
245 /*==========================================================================
246  = =
247  = Public methods =
248  = =
249  ==========================================================================*/
250 
251 /***********************************************************************//**
252  * @brief Returns MC instrument direction
253  *
254  * @param[in] energy Event energy.
255  * @param[in] time Event time.
256  * @param[in] obs CTA Observation.
257  * @param[in,out] ran Random number generator.
258  * @return Instrument direction
259  *
260  * @exception GException::invalid_return_value
261  * Spatial background model is empty, or method repeatedly
262  * encounters zero model values.
263  *
264  * Return random instrument direction.
265  ***************************************************************************/
267  const GTime& time,
268  const GCTAObservation& obs,
269  GRan& ran) const
270 {
271  // Initialise instrument direction
272  GCTAInstDir dir;
273 
274  // Get maximum value for Monte Carlo simulations (include some margin)
275  double mc_max_value = 1.5 * this->mc_max_value(obs);
276 
277  // Throw an exception if maximum value is zero
278  if (mc_max_value <= 0.0) {
279  std::string msg = "Spatial background model is empty. Please provide "
280  "a valid spatial background model.";
282  }
283 
284  // Get DETX and DETY value of RoI centre in radians
285  GCTARoi roi = obs.roi();
286  GCTAInstDir roi_centre = roi.centre();
287  double detx_centre = roi_centre.detx();
288  double dety_centre = roi_centre.dety();
289 
290  // Get DETX and DETY minima and maximum in radians
291  double radius = obs.roi().radius() * gammalib::deg2rad;
292  double detx_min = detx_centre - radius;
293  double detx_max = detx_centre + radius;
294  double dety_min = dety_centre - radius;
295  double dety_max = dety_centre + radius;
296  double detx_width = detx_max - detx_min;
297  double dety_width = dety_max - dety_min;
298 
299  // Initialise rejection method
300  double f = 0.0;
301  double ftest = 1.0;
302  int zeros = 0;
303  const int max_subsequent_zeros = 10;
304 
305  // Find instrument direction by rejection method
306  while (ftest > f) {
307 
308  // Draw random DETX/DETY until a non-zero function value is found.
309  // Subsequent zero values may indicate that there is an issue with
310  // the model. After a limited number of subsequent zeros the loop
311  // is terminated with an exception
312  int zeros = 0;
313  do {
314 
315  // Draw random DETX and DETY coordinates
316  double detx = detx_min + detx_width * ran.uniform();
317  double dety = dety_min + dety_width * ran.uniform();
318 
319  // Set instrument direction
320  dir.detx(detx);
321  dir.dety(dety);
322 
323  // Derive sky direction from instrument coordinates
324  GSkyDir skydir = obs.pointing().skydir(dir);
325 
326  // Set sky direction in GCTAInstDir object
327  dir.dir(skydir);
328 
329  // Compute function value if the coordinate is comprised in
330  // the RoI
331  if (roi.contains(dir)) {
332  f = eval(dir, energy, time, false);
333  }
334  else {
335  f = 0.0;
336  }
337 
338  // Handle zero function value
339  if (f == 0.0) {
340  zeros++;
341  if (zeros > max_subsequent_zeros) {
342  std::string msg = "Method repeatedly encounters zero values "
343  "of spatial background model. Make sure "
344  "that the background model is correct.";
346  }
347  }
348 
349  } while (f == 0);
350 
351  // Get uniform random value between zero and the maximum model value.
352  // The loop is quit if the random number is smaller than the function
353  // value.
354  ftest = ran.uniform() * mc_max_value;
355 
356  } // endwhile
357 
358  // Return instrument direction
359  return dir;
360 }
361 
362 
363 /***********************************************************************//**
364  * @brief Return integral of spatial model component
365  *
366  * @param[in] energy Measured event energy.
367  * @param[in] time Measured event time.
368  * @param[in] obs Observation.
369  * @return Integral of spatial model component.
370  *
371  * Spatially integrates the spatial background model component for a given
372  * measured event energy and event time over the region of interest (RoI)
373  * of a given observation.
374  *
375  * The method uses a 2D Romberg integration to numerically integrate the
376  * spatial background model component.
377  ***************************************************************************/
378 double GCTAModelSpatial::npred(const GEnergy& energy,
379  const GTime& time,
380  const GObservation& obs) const
381 {
382  // Set number of iterations for Romberg integration.
383  static const int min_iter_theta = 5;
384  static const int min_iter_phi = 5;
385  static const int max_iter_theta = 8;
386  static const int max_iter_phi = 8;
387 
388  // Initialise result
389  double npred = 0.0;
390 
391  // Get reference on CTA pointing and event list from observation
392  const GCTAPointing& pnt = gammalib::cta_pnt(G_NPRED, obs);
393  const GCTAEventList& events = gammalib::cta_event_list(G_NPRED, obs);
394 
395  // Get instrument direction of RoI centre
396  GCTAInstDir roi_centre = pnt.instdir(events.roi().centre().dir());
397 
398  // Get RoI radius in radians
399  double roi_radius = events.roi().radius() * gammalib::deg2rad;
400 
401  // Set number of radial integration iterations
402  int iter_theta = gammalib::iter_rho(roi_radius, g_npred_resolution,
403  min_iter_theta, max_iter_theta);
404 
405  // Setup integration function
407  energy,
408  time,
409  roi_centre,
410  min_iter_phi,
411  max_iter_phi);
412 
413  // Setup integration
414  GIntegral integral(&integrand);
415 
416  // Set fixed number of iterations
417  integral.fixed_iter(iter_theta);
418 
419  // Spatially integrate radial component (assumes that RoI centre is
420  // at DETX=DETY=0)
421  npred = integral.romberg(0.0, roi_radius);
422 
423  // Debug: Check for NaN
424  #if defined(G_NAN_CHECK)
425  if (gammalib::is_notanumber(npred) || gammalib::is_infinite(npred)) {
426  std::string origin = "GCTAModelSpatial::npred";
427  std::string message = " NaN/Inf encountered (npred=" +
428  gammalib::str(npred) + ", roi_radius=" +
429  gammalib::str(roi_radius) + ")";
430  gammalib::warning(origin, message);
431  }
432  #endif
433 
434  // Return Npred
435  return npred;
436 }
437 
438 
439 /*==========================================================================
440  = =
441  = Private methods =
442  = =
443  ==========================================================================*/
444 
445 /***********************************************************************//**
446  * @brief Initialise class members
447  ***************************************************************************/
449 {
450  // Initialise members
451  m_pars.clear();
452 
453  // Return
454  return;
455 }
456 
457 
458 /***********************************************************************//**
459  * @brief Copy class members
460  *
461  * @param[in] model Spatial background model component.
462  ***************************************************************************/
464 {
465  // Copy members
466  m_pars = model.m_pars;
467 
468  // Return
469  return;
470 }
471 
472 
473 /***********************************************************************//**
474  * @brief Delete class members
475  ***************************************************************************/
477 {
478  // Return
479  return;
480 }
481 
482 
483 /***********************************************************************//**
484  * @brief Kernel for offset angle integration of spatial component
485  *
486  * @param[in] theta Offset angle from RoI centre (radians).
487  *
488  * Computes
489  *
490  * \f[
491  * K(\rho | E, t) = \sin \theta \times
492  * \int_{0}^{2\pi}
493  * B(\theta,\phi | E, t) d\phi
494  * \f]
495  *
496  * where \f$B(\theta,\phi | E, t)\f$ is the spatial component of the
497  * background model for a specific observed energy \f$E\f$ and time \f$t\f$.
498  ***************************************************************************/
500 {
501  // Initialise value
502  double value = 0.0;
503 
504  // Continue only if offset angle is positive
505  if (theta > 0.0) {
506 
507  // Setup phi integration kernel
509  m_energy,
510  m_time,
511  m_roi_centre,
512  theta);
513 
514  // Setup integration
515  GIntegral integral(&integrand);
516 
517  // Set number of azimuthal integration iterations
518  int iter = gammalib::iter_phi(theta, g_npred_resolution,
520 
521  // Set fixed number of iterations
522  integral.fixed_iter(iter);
523 
524  // Integrate over phi
525  value = integral.romberg(0.0, gammalib::twopi) * std::sin(theta);
526 
527  // Debug: Check for NaN
528  #if defined(G_NAN_CHECK)
529  if (gammalib::is_notanumber(value) || gammalib::is_infinite(value)) {
530  std::string origin = "GCTAModelSpatial::npred_roi_kern_theta::eval"
531  "(" + gammalib::str(theta) + ")";
532  std::string message = " NaN/Inf encountered (value=" +
533  gammalib::str(value) + ")";
534  gammalib::warning(origin, message);
535  }
536  #endif
537 
538  } // endif: offset angle was positive
539 
540  // Return value
541  return value;
542 }
543 
544 
545 /***********************************************************************//**
546  * @brief Kernel for azimuth angle integration of spatial component
547  *
548  * @param[in] phi Azimuth angle around RoI centre (radians).
549  *
550  * Computes
551  *
552  * \f[
553  * B(\theta, \phi | E, t)
554  * \f]
555  *
556  * using
557  *
558  * \f[ {\rm detx} = \theta \cos \phi \f]
559  * \f[ {\rm dety} = \theta \sin \phi \f]
560  *
561  * @todo Verify correct orientation of detx and dety with respect to phi
562  ***************************************************************************/
564 {
565  // Compute detx and dety in radians
566  double detx = m_roi_centre.detx();
567  double dety = m_roi_centre.dety();
568  if (m_theta > 0.0 ) {
569  detx += m_theta * std::cos(phi);
570  dety += m_theta * std::sin(phi);
571  }
572 
573  // Setup CTA instrument direction
574  GCTAInstDir dir(detx, dety);
575 
576  // Get background value
577  double value = m_spatial->eval(dir, m_energy, m_time);
578 
579  // Debug: Check for NaN
580  #if defined(G_NAN_CHECK)
581  if (gammalib::is_notanumber(value) || gammalib::is_infinite(value)) {
582  std::string origin = "GCTAModelSpatial::npred_roi_kern_phi::eval"
583  "(" + gammalib::str(phi) + ")";
584  std::string message = " NaN/Inf encountered (value=" +
585  gammalib::str(value) + ", detx=" +
586  gammalib::str(detx) + ", dety=" +
587  gammalib::str(dety) + ")";
588  gammalib::warning(origin, message);
589  }
590  #endif
591 
592  // Return Npred
593  return value;
594 }
virtual GCTAInstDir mc(const GEnergy &energy, const GTime &time, const GCTAObservation &obs, GRan &ran) const
Returns MC instrument direction.
Abstract spatial model class.
CTA response helper classes definition.
#define G_NPRED
CTA event list class interface definition.
double romberg(std::vector< double > bounds, const int &order=5)
Perform Romberg integration.
Definition: GIntegral.cpp:380
void detx(const double &x)
Set DETX coordinate (in radians)
void warning(const std::string &origin, const std::string &message)
Emits warning.
Definition: GTools.cpp:1265
const double g_npred_resolution
Scale of bkg.
virtual void roi(const GRoi &roi)
Set ROI.
void copy_members(const GCTAModelSpatial &model)
Copy class members.
double eval(const double &phi)
Kernel for azimuth angle integration of spatial component.
GVector cos(const GVector &vector)
Computes cosine of vector elements.
Definition: GVector.cpp:1100
GCTAModelSpatial(void)
Void constructor.
const GCTAPointing & cta_pnt(const std::string &origin, const GObservation &obs)
Retrieve CTA pointing from generic observation.
virtual double mc_max_value(const GCTAObservation &obs) const =0
CTA event list class.
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
Random number generator class.
Definition: GRan.hpp:44
Interface for the CTA region of interest class.
Definition: GCTARoi.hpp:48
Time class.
Definition: GTime.hpp:54
Gammalib tools definition.
GIntegral class interface definition.
Definition: GIntegral.hpp:45
void pointing(const GCTAPointing &pointing)
Set CTA pointing.
void dir(const GSkyDir &dir)
Set sky direction.
int m_min_iter
Minimum number of Romberg iterations.
bool is_notanumber(const double &x)
Signal if argument is not a number.
Definition: GTools.hpp:184
bool is_infinite(const double &x)
Signal if argument is infinite.
Definition: GTools.hpp:167
Definition of support function used by CTA classes.
Model parameter class.
Definition: GModelPar.hpp:87
int m_max_iter
Maximum number of Romberg iterations.
const double deg2rad
Definition: GMath.hpp:43
#define G_ACCESS1
virtual ~GCTAModelSpatial(void)
Destructor.
GCTARoi roi(void) const
Get Region of Interest.
#define G_ACCESS2
const GCTAInstDir & centre(void) const
Returns region of interest centre.
Definition: GCTARoi.hpp:110
void free_members(void)
Delete class members.
const GCTAModelSpatial * m_spatial
Pointer to spatial component.
GCTAInstDir instdir(const GSkyDir &skydir) const
Get instrument direction from sky direction.
double uniform(void)
Returns random double precision floating value in range 0 to 1.
Definition: GRan.cpp:242
void init_members(void)
Initialise class members.
CTA pointing class.
CTA instrument direction class interface definition.
const GCTAEventList & cta_event_list(const std::string &origin, const GObservation &obs)
Retrieve CTA event list from generic observation.
virtual bool contains(const GEvent &event) const
Check if region of interest contains an event.
Definition: GCTARoi.cpp:214
int size(void) const
Return number of model parameters.
virtual GModelPar & operator[](const int &index)
Returns model parameter.
Abstract observation base class.
Abstract observation base class interface definition.
int iter_phi(const double &rho, const double &resolution, const int &iter_min, const int &iter_max)
Determine number of azimuthal Romberg iterations.
Abstract spatial model class interface definition.
GVector sin(const GVector &vector)
Computes sine of vector elements.
Definition: GVector.cpp:1226
Exception handler interface definition.
virtual double eval(const GCTAInstDir &dir, const GEnergy &energy, const GTime &time, const bool &gradients=false) const =0
void fixed_iter(const int &iter)
Set fixed number of iterations.
Definition: GIntegral.hpp:181
#define G_MC
CTA instrument direction class.
Definition: GCTAInstDir.hpp:59
void dety(const double &y)
Set DETY coordinate (in radians)
const double twopi
Definition: GMath.hpp:36
CTA observation class.
std::vector< GModelPar * > m_pars
Parameter pointers.
virtual double npred(const GEnergy &energy, const GTime &time, const GObservation &obs) const
Return integral of spatial model component.
double eval(const double &theta)
Kernel for offset angle integration of spatial component.
Integration class interface definition.
int iter_rho(const double &rho_max, const double &resolution, const int &iter_min, const int &iter_max)
Determine number of radial Romberg iterations.
Sky direction class.
Definition: GSkyDir.hpp:62
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:413