GammaLib  1.7.0.dev
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
GModelSpatialPointSource.cpp
Go to the documentation of this file.
1 /***************************************************************************
2  * GModelSpatialPointSource.cpp - Spatial point source model class *
3  * ----------------------------------------------------------------------- *
4  * copyright (C) 2009-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 GModelSpatialPointSource.cpp
23  * @brief Point source 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 "GMath.hpp"
36 
37 /* __ Constants __________________________________________________________ */
38 const double tolerance = 0.000027777778; // angular tolerance is 1 arcsec
39 
40 
41 /* __ Globals ____________________________________________________________ */
43 const GModelSpatialRegistry g_spatial_ptsrc_registry(&g_spatial_ptsrc_seed);
44 #if defined(G_LEGACY_XML_FORMAT)
45 const GModelSpatialPointSource g_spatial_ptsrc_legacy_seed(true, "SkyDirFunction");
46 const GModelSpatialRegistry g_spatial_ptsrc_legacy_registry(&g_spatial_ptsrc_legacy_seed);
47 #endif
48 
49 /* __ Method name definitions ____________________________________________ */
50 #define G_READ "GModelSpatialPointSource::read(GXmlElement&)"
51 #define G_WRITE "GModelSpatialPointSource::write(GXmlElement&)"
52 
53 /* __ Macros _____________________________________________________________ */
54 
55 /* __ Coding definitions _________________________________________________ */
56 
57 /* __ Debug definitions __________________________________________________ */
58 
59 
60 /*==========================================================================
61  = =
62  = Constructors/destructors =
63  = =
64  ==========================================================================*/
65 
66 /***********************************************************************//**
67  * @brief Void constructor
68  *
69  * Constructs empty point source model.
70  ***************************************************************************/
72 {
73  // Initialise members
74  init_members();
75 
76  // Return
77  return;
78 }
79 
80 
81 /***********************************************************************//**
82  * @brief Model type constructor
83  *
84  * @param[in] dummy Dummy flag.
85  * @param[in] type Model type.
86  *
87  * Constructs empty point source model by specifying a model @p type.
88  ***************************************************************************/
90  const std::string& type) :
92 {
93  // Initialise members
94  init_members();
95 
96  // Set model type
97  m_type = type;
98 
99  // Return
100  return;
101 }
102 
103 
104 /***********************************************************************//**
105  * @brief Sky direction constructor
106  *
107  * @param[in] dir Sky direction.
108  *
109  * Construct a point source spatial model from a sky direction.
110  ***************************************************************************/
112  GModelSpatial()
113 {
114  // Initialise members
115  init_members();
116 
117  // Assign direction
118  this->dir(dir);
119 
120  // Return
121  return;
122 }
123 
124 
125 /***********************************************************************//**
126  * @brief Value constructor
127  *
128  * @param[in] ra Right Ascencion of model centre.
129  * @param[in] dec Declination of model centre.
130  *
131  * Construct a point source spatial model from the Right Ascension and
132  * Declination of the model centre.
133  ***************************************************************************/
135  const double& dec) :
136  GModelSpatial()
137 {
138  // Initialise members
139  init_members();
140 
141  // Set values
142  m_ra.value(ra);
143  m_dec.value(dec);
144 
145  // Return
146  return;
147 }
148 
149 
150 /***********************************************************************//**
151  * @brief XML constructor
152  *
153  * @param[in] xml XML element.
154  *
155  * Construct a point source spatial model by extracting information from an
156  * XML element. See the read() method for more information about the expected
157  * structure of the XML element.
158  ***************************************************************************/
160  GModelSpatial()
161 {
162  // Initialise members
163  init_members();
164 
165  // Read information from XML element
166  read(xml);
167 
168  // Return
169  return;
170 }
171 
172 
173 /***********************************************************************//**
174  * @brief Copy constructor
175  *
176  * @param[in] model Point source spatial model.
177  ***************************************************************************/
179  GModelSpatial(model)
180 {
181  // Initialise members
182  init_members();
183 
184  // Copy members
185  copy_members(model);
186 
187  // Return
188  return;
189 }
190 
191 
192 /***********************************************************************//**
193  * @brief Destructor
194  ***************************************************************************/
196 {
197  // Free members
198  free_members();
199 
200  // Return
201  return;
202 }
203 
204 
205 /*==========================================================================
206  = =
207  = Operators =
208  = =
209  ==========================================================================*/
210 
211 /***********************************************************************//**
212  * @brief Assignment operator
213  *
214  * @param[in] model Point source spatial model.
215  * @return Point source spatial model.
216  ***************************************************************************/
218 {
219  // Execute only if object is not identical
220  if (this != &model) {
221 
222  // Copy base class members
223  this->GModelSpatial::operator=(model);
224 
225  // Free members
226  free_members();
227 
228  // Initialise members
229  init_members();
230 
231  // Copy members
232  copy_members(model);
233 
234  } // endif: object was not identical
235 
236  // Return
237  return *this;
238 }
239 
240 
241 /*==========================================================================
242  = =
243  = Public methods =
244  = =
245  ==========================================================================*/
246 
247 /***********************************************************************//**
248  * @brief Clear point source model
249  ***************************************************************************/
251 {
252  // Free class members (base and derived classes, derived class first)
253  free_members();
255 
256  // Initialise members
258  init_members();
259 
260  // Return
261  return;
262 }
263 
264 
265 /***********************************************************************//**
266  * @brief Clone point source model
267  *
268  * @return Pointer to deep copy of point source model.
269  ***************************************************************************/
271 {
272  // Clone point source model
273  return new GModelSpatialPointSource(*this);
274 }
275 
276 
277 /***********************************************************************//**
278  * @brief Evaluate function
279  *
280  * @param[in] photon Incident photon.
281  * @param[in] gradients Compute gradients?
282  * @return Model value.
283  *
284  * Evaluates the spatial part for a point source model. It implements a delta
285  * function with respect to the coordinates of the source. For numerical
286  * reasons, a certain tolerance is accepted (typically 0.1 arcsec, i.e. well
287  * below the angular resolution of gamma-ray telescopes).
288  *
289  * The method will not compute analytical parameter gradients, even if the
290  * @p gradients argument is set to true. Point source parameter gradients
291  * need to be computed numerically.
292  ***************************************************************************/
294  const bool& gradients) const
295 {
296  // Set value dependent on source distance
297  double value = (photon.dir().dist_deg(dir()) < tolerance) ? 1.0 : 0.0;
298 
299  // Return value
300  return value;
301 }
302 
303 
304 /***********************************************************************//**
305  * @brief Returns MC sky direction
306  *
307  * @param[in] energy Photon energy.
308  * @param[in] time Photon arrival time.
309  * @param[in,out] ran Random number generator.
310  * @return Sky direction.
311  *
312  * Draws an arbitrary sky direction for the point source model. As the point
313  * source is a point in the sky, the methods always returns the directon of
314  * the point source.
315  ***************************************************************************/
317  const GTime& time,
318  GRan& ran) const
319 {
320  // Return sky direction
321  return (dir());
322 }
323 
324 
325 /***********************************************************************//**
326  * @brief Checks where model contains specified sky direction
327  *
328  * @param[in] dir Sky direction.
329  * @param[in] margin Margin to be added to sky direction (degrees)
330  * @return True if the model contains the sky direction.
331  *
332  * Signals whether a sky direction is contained in the point source
333  * model.
334  ***************************************************************************/
336  const double& margin) const
337 {
338  // Compute distance to centre (radian)
339  double distance = dir.dist(this->dir());
340 
341  // Return flag
342  return (distance <= margin*gammalib::deg2rad);
343 }
344 
345 
346 /***********************************************************************//**
347  * @brief Read model from XML element
348  *
349  * @param[in] xml XML element containing point source model information.
350  *
351  * @exception GException::model_invalid_parnum
352  * Invalid number of model parameters found in XML element.
353  * @exception GException::model_invalid_parnames
354  * Invalid model parameter names found in XML element.
355  *
356  * Read the point source information from an XML element with the following
357  * format
358  *
359  * <spatialModel type="PointSource">
360  * <parameter free="0" max="360" min="-360" name="RA" scale="1" value="83.6331" />
361  * <parameter free="0" max="90" min="-90" name="DEC" scale="1" value="22.0145" />
362  * </spatialModel>
363  *
364  * or
365  *
366  * <spatialModel type="PointSource">
367  * <parameter free="0" max="360" min="-360" name="GLON" scale="1" value="83.6331" />
368  * <parameter free="0" max="90" min="-90" name="GLAT" scale="1" value="22.0145" />
369  * </spatialModel>
370  *
371  ***************************************************************************/
373 {
374  // Read RA/DEC parameters
375  if (gammalib::xml_has_par(xml, "RA") && gammalib::xml_has_par(xml, "DEC")) {
376 
377  // Get parameters
378  const GXmlElement* ra = gammalib::xml_get_par(G_READ, xml, "RA");
379  const GXmlElement* dec = gammalib::xml_get_par(G_READ, xml, "DEC");
380 
381  // Read parameters
382  m_ra.read(*ra);
383  m_dec.read(*dec);
384 
385  }
386 
387  // ... otherwise read GLON/GLAT parameters
388  else {
389 
390  // Get parameters
391  const GXmlElement* glon = gammalib::xml_get_par(G_READ, xml, "GLON");
392  const GXmlElement* glat = gammalib::xml_get_par(G_READ, xml, "GLAT");
393 
394  // Read parameters
395  m_ra.read(*glon);
396  m_dec.read(*glat);
397 
398  // Convert into RA/DEC
399  GSkyDir dir;
400  dir.lb_deg(ra(), dec()),
401  m_ra.value(dir.ra_deg());
402  m_dec.value(dir.dec_deg());
403 
404  // Set names to RA/DEC
405  m_ra.name("RA");
406  m_dec.name("DEC");
407 
408  }
409 
410  // Return
411  return;
412 }
413 
414 
415 /***********************************************************************//**
416  * @brief Write model into XML element
417  *
418  * @param[in] xml XML element into which model information is written.
419  *
420  * @exception GException::model_invalid_spatial
421  * Existing XML element is not of type 'SkyDirFunction'
422  * @exception GException::model_invalid_parnum
423  * Invalid number of model parameters found in XML element.
424  * @exception GException::model_invalid_parnames
425  * Invalid model parameter names found in XML element.
426  *
427  * Write the point source information into an XML element with the following
428  * format
429  *
430  * <spatialModel type="PointSource">
431  * <parameter free="0" max="360" min="-360" name="RA" scale="1" value="83.6331" />
432  * <parameter free="0" max="90" min="-90" name="DEC" scale="1" value="22.0145" />
433  * </spatialModel>
434  *
435  * @todo The case that an existing spatial XML element with "GLON" and "GLAT"
436  * as coordinates is not supported.
437  ***************************************************************************/
439 {
440  // Set model type
441  if (xml.attribute("type") == "") {
442  xml.attribute("type", type());
443  }
444 
445  // Verify model type
446  if (xml.attribute("type") != type()) {
448  "Spatial model is not of type \""+type()+"\".");
449  }
450 
451  // Get or create parameters
454 
455  // Write parameters
456  m_ra.write(*ra);
457  m_dec.write(*dec);
458 
459  // Return
460  return;
461 }
462 
463 
464 /***********************************************************************//**
465  * @brief Print point source information
466  *
467  * @param[in] chatter Chattiness.
468  * @return String containing point source information.
469  ***************************************************************************/
470 std::string GModelSpatialPointSource::print(const GChatter& chatter) const
471 {
472  // Initialise result string
473  std::string result;
474 
475  // Continue only if chatter is not silent
476  if (chatter != SILENT) {
477 
478  // Append header
479  result.append("=== GModelSpatialPointSource ===");
480 
481  // Append parameters
482  result.append("\n"+gammalib::parformat("Number of parameters"));
483  result.append(gammalib::str(size()));
484  for (int i = 0; i < size(); ++i) {
485  result.append("\n"+m_pars[i]->print(chatter));
486  }
487 
488  } // endif: chatter was not silent
489 
490  // Return result
491  return result;
492 }
493 
494 
495 /***********************************************************************//**
496  * @brief Return position of point source
497  *
498  * @return Point source sky direction.
499  *
500  * Returns the sky direction of the point source.
501  ***************************************************************************/
503 {
504  // Allocate sky direction
505  GSkyDir srcDir;
506 
507  // Set sky direction
508  srcDir.radec_deg(ra(), dec());
509 
510  // Return direction
511  return srcDir;
512 }
513 
514 
515 /***********************************************************************//**
516  * @brief Set position of point source
517  *
518  * Sets the sky direction of the point source.
519  ***************************************************************************/
521 {
522  // Assign Right Ascension and Declination
523  m_ra.value(dir.ra_deg());
524  m_dec.value(dir.dec_deg());
525 
526  // Return
527  return;
528 }
529 
530 
531 /*==========================================================================
532  = =
533  = Private methods =
534  = =
535  ==========================================================================*/
536 
537 /***********************************************************************//**
538  * @brief Initialise class members
539  ***************************************************************************/
541 {
542  // Initialise model type
543  m_type = "PointSource";
544 
545  // Initialise Right Ascension
546  m_ra.clear();
547  m_ra.name("RA");
548  m_ra.unit("deg");
549  m_ra.fix();
550  m_ra.scale(1.0);
551  m_ra.gradient(0.0);
552  m_ra.has_grad(false);
553 
554  // Initialise Declination
555  m_dec.clear();
556  m_dec.name("DEC");
557  m_dec.unit("deg");
558  m_dec.fix();
559  m_dec.scale(1.0);
560  m_dec.gradient(0.0);
561  m_dec.has_grad(false);
562 
563  // Set parameter pointer(s)
564  m_pars.clear();
565  m_pars.push_back(&m_ra);
566  m_pars.push_back(&m_dec);
567 
568  // Initialise other members
569  m_region.clear();
570 
571  // Return
572  return;
573 }
574 
575 
576 /***********************************************************************//**
577  * @brief Copy class members
578  *
579  * @param[in] model Point source spatial model.
580  ***************************************************************************/
582 {
583  // Copy members
584  m_type = model.m_type;
585  m_ra = model.m_ra;
586  m_dec = model.m_dec;
587  m_region = model.m_region;
588 
589  // Set parameter pointer(s)
590  m_pars.clear();
591  m_pars.push_back(&m_ra);
592  m_pars.push_back(&m_dec);
593 
594  // Return
595  return;
596 }
597 
598 
599 /***********************************************************************//**
600  * @brief Delete class members
601  ***************************************************************************/
603 {
604  // Return
605  return;
606 }
607 
608 
609 /***********************************************************************//**
610  * @brief Set boundary sky region
611  ***************************************************************************/
613 {
614  // Set sky region centre
616 
617  // Set sky region radius to zero
618  m_region.radius(0.0);
619 
620  // Return
621  return;
622 }
void free_members(void)
Delete class members.
double dist_deg(const GSkyDir &dir) const
Compute angular distance between sky directions in degrees.
Definition: GSkyDir.hpp:280
double dec_deg(void) const
Returns Declination in degrees.
Definition: GSkyDir.hpp:250
Point source spatial model class interface definition.
const std::string & name(void) const
Return parameter name.
double gradient(void) const
Return parameter gradient.
int size(void) const
Return number of parameters.
void write(GXmlElement &xml) const
Set or update parameter attributes in XML element.
Definition: GModelPar.cpp:353
double dec(void) const
Return Declination of model centre.
XML element node class.
Definition: GXmlElement.hpp:47
virtual GModelSpatialPointSource * clone(void) const
Clone point source model.
virtual GModelSpatial & operator=(const GModelSpatial &model)
Assignment operator.
virtual double eval(const GPhoton &photon, const bool &gradients=false) const
Evaluate function.
Random number generator class.
Definition: GRan.hpp:44
Time class.
Definition: GTime.hpp:54
Gammalib tools definition.
double ra(void) const
Return Right Ascencion of model centre.
virtual void clear(void)
Clear point source model.
Interface definition for the spatial model registry class.
void copy_members(const GModelSpatialPointSource &model)
Copy class members.
const double & radius(void) const
Return circular region radius (in degrees)
Class that handles photons.
Definition: GPhoton.hpp:47
const double & scale(void) const
Return parameter scale.
const double deg2rad
Definition: GMath.hpp:43
bool has_grad(void) const
Signal if parameter gradient is computed analytically.
virtual std::string print(const GChatter &chatter=NORMAL) const
Print point source information.
const GXmlAttribute * attribute(const int &index) const
Return attribute.
virtual ~GModelSpatialPointSource(void)
Destructor.
virtual void write(GXmlElement &xml) const
Write model into XML element.
virtual GSkyDir mc(const GEnergy &energy, const GTime &time, GRan &ran) const
Returns MC sky direction.
std::vector< GModelPar * > m_pars
Parameter pointers.
const double tolerance
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:1513
virtual void read(const GXmlElement &xml)
Read model from XML element.
void clear(void)
Clear parameter.
Spatial model registry class definition.
double ra_deg(void) const
Returns Right Ascension in degrees.
Definition: GSkyDir.hpp:223
GChatter
Definition: GTypemaps.hpp:33
void radec_deg(const double &ra, const double &dec)
Set equatorial sky direction (degrees)
Definition: GSkyDir.cpp:202
const GModelSpatialPointSource g_spatial_ptsrc_seed
Point source spatial model.
bool xml_has_par(const GXmlElement &xml, const std::string &name)
Checks if parameter with given name in XML element exists.
Definition: GTools.cpp:1475
virtual std::string type(void) const
Return model type.
const GSkyDir & centre(void) const
Return circular region centre.
void free_members(void)
Delete class members.
GSkyDir dir(void) const
Return position of point source.
GModelPar m_ra
Right Ascension (deg)
virtual bool contains(const GSkyDir &dir, const double &margin=0.0) const
Checks where model contains specified sky direction.
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.
void lb_deg(const double &l, const double &b)
Set galactic sky direction (degrees)
Definition: GSkyDir.cpp:256
Exception handler interface definition.
double dist(const GSkyDir &dir) const
Compute angular distance between sky directions in radians.
Definition: GSkyDir.hpp:265
GModelSpatialPointSource(void)
Void constructor.
void init_members(void)
Initialise class members.
Abstract spatial model base class.
GModelPar m_dec
Declination (deg)
void set_region(void) const
Set boundary sky region.
#define G_WRITE
void init_members(void)
Initialise class members.
virtual GModelSpatialPointSource & operator=(const GModelSpatialPointSource &model)
Assignment operator.
std::string parformat(const std::string &s, const int &indent=0)
Convert string in parameter format.
Definition: GTools.cpp:1022
#define G_READ
GSkyRegionCircle m_region
Bounding circle.
Sky direction class.
Definition: GSkyDir.hpp:62
const GSkyDir & dir(void) const
Return photon sky direction.
Definition: GPhoton.hpp:110
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:1562
Mathematical function definitions.
void clear(void)
Clear instance.
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