GammaLib  1.7.0.dev
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
GModelSpatialDiffuseConst.cpp
Go to the documentation of this file.
1 /***************************************************************************
2  * GModelSpatialDiffuseConst.cpp - Spatial isotropic model class *
3  * ----------------------------------------------------------------------- *
4  * copyright (C) 2010-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 GModelSpatialDiffuseConst.cpp
23  * @brief Isotropic 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 
39 /* __ Globals ____________________________________________________________ */
41 const GModelSpatialRegistry g_spatial_const_registry(&g_spatial_const_seed);
42 #if defined(G_LEGACY_XML_FORMAT)
43 const GModelSpatialDiffuseConst g_spatial_const_legacy_seed(true, "ConstantValue");
44 const GModelSpatialRegistry g_spatial_const_legacy_registry(&g_spatial_const_legacy_seed);
45 #endif
46 
47 /* __ Method name definitions ____________________________________________ */
48 #define G_MC_NORM "GModelSpatialDiffuseConst::mc_norm(GSkyDir&, double&)"
49 #define G_READ "GModelSpatialDiffuseConst::read(GXmlElement&)"
50 #define G_WRITE "GModelSpatialDiffuseConst::write(GXmlElement&)"
51 
52 /* __ Macros _____________________________________________________________ */
53 
54 /* __ Coding definitions _________________________________________________ */
55 
56 /* __ Debug definitions __________________________________________________ */
57 
58 
59 /*==========================================================================
60  = =
61  = Constructors/destructors =
62  = =
63  ==========================================================================*/
64 
65 /***********************************************************************//**
66  * @brief Void constructor
67  *
68  * Constructs empty isotropic spatial model.
69  ***************************************************************************/
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 isotropic spatial 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 XML constructor
106  *
107  * @param[in] xml XML element.
108  *
109  * Constructs isotropic spatial model by extracting information from an XML
110  * element. See the read() method for more information about the expected
111  * structure of the XML element.
112  ***************************************************************************/
115 {
116  // Initialise members
117  init_members();
118 
119  // Read information from XML element
120  read(xml);
121 
122  // Return
123  return;
124 }
125 
126 
127 /***********************************************************************//**
128  * @brief Value constructor
129  *
130  * @param[in] value Isotropic intensity value (\f$sr^{-1}\f$).
131  *
132  * Constructs isotropic spatial model by assigning an intensity of @p value
133  * in units of \f$sr^{-1}\f$ to the diffuse emission. This constructor
134  * explicitly sets the @a m_value member of the model.
135  ***************************************************************************/
138 {
139  // Initialise members
140  init_members();
141 
142  // Set parameter
143  m_value.value(value);
144 
145  // Perform autoscaling of parameter
146  autoscale();
147 
148  // Return
149  return;
150 }
151 
152 
153 /***********************************************************************//**
154  * @brief Copy constructor
155  *
156  * @param[in] model Isotropic spatial model.
157  *
158  * Constructs isotropic spatial model by copying another isotropic spatial
159  * model.
160  ***************************************************************************/
162  GModelSpatialDiffuse(model)
163 {
164  // Initialise members
165  init_members();
166 
167  // Copy members
168  copy_members(model);
169 
170  // Return
171  return;
172 }
173 
174 
175 /***********************************************************************//**
176  * @brief Destructor
177  *
178  * Destructs isotropic spatial model.
179  ***************************************************************************/
181 {
182  // Free members
183  free_members();
184 
185  // Return
186  return;
187 }
188 
189 
190 /*==========================================================================
191  = =
192  = Operators =
193  = =
194  ==========================================================================*/
195 
196 /***********************************************************************//**
197  * @brief Assignment operator
198  *
199  * @param[in] model Isotropic spatial model.
200  * @return Isotropic spatial model.
201  *
202  * Assigns an isotropic spatial model.
203  ***************************************************************************/
205 {
206  // Execute only if object is not identical
207  if (this != &model) {
208 
209  // Copy base class members
210  this->GModelSpatialDiffuse::operator=(model);
211 
212  // Free members
213  free_members();
214 
215  // Initialise members
216  init_members();
217 
218  // Copy members
219  copy_members(model);
220 
221  } // endif: object was not identical
222 
223  // Return
224  return *this;
225 }
226 
227 
228 /*==========================================================================
229  = =
230  = Public methods =
231  = =
232  ==========================================================================*/
233 
234 /***********************************************************************//**
235  * @brief Clear isotropic spatial model
236  *
237  * Clears the isotropic spatial model. This method is equivalent to creating
238  * an isotropic spatial model using the void constructor.
239  ***************************************************************************/
241 {
242  // Free class members
243  free_members();
246 
247  // Initialise members
250  init_members();
251 
252  // Return
253  return;
254 }
255 
256 
257 /***********************************************************************//**
258  * @brief Clone isotropic spatial model
259  *
260  * @return Pointer to deep copy of isotropic spatial model.
261  *
262  * Returns a pointer to a deep copy of the isotropic spatial model.
263  ***************************************************************************/
265 {
266  // Clone isotropic spatial model
267  return new GModelSpatialDiffuseConst(*this);
268 }
269 
270 
271 /***********************************************************************//**
272  * @brief Evaluate isotropic spatial model value
273  *
274  * @param[in] photon Incident photon.
275  * @param[in] gradients Compute gradients?
276  * @return Model value (\f$sr^{-1}\f$).
277  *
278  * Evaluates the isotropic spatial model for a given @p photon, characterised
279  * by a sky direction, energy and arrival time \f$(\vec{p},E,t)\f$. By
280  * definition, the isotropic spatial model and gradient are independent of
281  * sky direction, energy and time. The model value is given by
282  *
283  * \f[
284  * M_\mathrm{S}(\vec{p}|E,t) = \frac{v}{4\pi}
285  * \f]
286  *
287  * where \f$v\f$ is the value() parameter, which is divided by the solid
288  * angle \f$4\pi\f$ of the celestial sphere.
289  *
290  * If the @p gradients flag is true the method will also evaluate the partial
291  * derivatives of the model. The value() gradient is given by
292  *
293  * \f[
294  * \frac{\partial M_\mathrm{S}(\vec{p}|E,t)}{\partial v_v} =
295  * \frac{v_s}{4\pi}
296  * \f]
297  *
298  * where \f$v=v_v v_s\f$ is the factorisation of the value() parameter into
299  * a value \f$v_v\f$ and scale \f$v_s\f$ term.
300  ***************************************************************************/
302  const bool& gradients) const
303 {
304  // Set normalization constant
305  const double norm = 1.0 / gammalib::fourpi;
306 
307  // Compute model value
308  double value = norm * m_value.value();
309 
310  // Optionally compute partial derivatives
311  if (gradients) {
312 
313  // Compute partial derivatives of the parameter values
314  double g_norm = (m_value.is_free()) ? norm * m_value.scale() : 0.0;
315 
316  // Set gradient
317  m_value.factor_gradient(g_norm);
318 
319  } // endif: computed partial derivatives
320 
321  // Return model value
322  return (value);
323 }
324 
325 
326 /***********************************************************************//**
327  * @brief Return MC sky direction
328  *
329  * @param[in] energy Photon energy.
330  * @param[in] time Photon arrival time.
331  * @param[in,out] ran Random number generator.
332  * @return Sky direction.
333  *
334  * Returns an arbitrary direction on the celestial sphere within a simulation
335  * cone region that has been defined by a previous call to the mc_norm()
336  * method (if no such call was made the entire sky will be assumed as the
337  * simulation cone). The sky direction is independent of event @p energy and
338  * arrival @p time.
339  ***************************************************************************/
341  const GTime& time,
342  GRan& ran) const
343 {
344  // Simulate offset from simulation cone centre
345  double theta = std::acos(1.0 - ran.uniform() * (1.0 - m_mc_cos_radius)) *
347  double phi = 360.0 * ran.uniform();
348 
349  // Initialise sky direction with simulation cone centre
350  GSkyDir dir(m_mc_centre);
351 
352  // Rotate sky direction by offset angles (phi, theta)
353  dir.rotate_deg(phi, theta);
354 
355  // Return sky direction
356  return dir;
357 }
358 
359 
360 /***********************************************************************//**
361  * @brief Return normalization of constant diffuse source for Monte Carlo
362  * simulations
363  *
364  * @param[in] dir Centre of simulation cone.
365  * @param[in] radius Radius of simulation cone (deg).
366  * @return Model normalization.
367  *
368  * @exception GException::invalid_argument
369  * Radius of simulation cone not coprised in [0,180] degrees.
370  *
371  * Returns the normalization of a constant diffuse source. The normalization
372  * is given by the product between the fraction \f$f\f$ of the sky covered
373  * by the simulation cone
374  *
375  * \f[
376  * f = \frac{2\pi \left( 1-\cos(r) \right)}{4\pi}
377  * \f]
378  *
379  * (where \f$r\f$ is the radius of the simulation cone) and the model
380  * parameter value(). The normalization only depends of the radius of the
381  * simulation cone and is invariant to its centre.
382  ***************************************************************************/
384  const double& radius) const
385 {
386  // Throw exception if the radius is invalid
387  if (radius < 0.0 || radius > 180.0) {
388  std::string msg = "Invalid simulation cone radius "+
389  gammalib::str(radius)+" specified. Please provide "
390  "a value comprised within 0 and 180 degrees.";
392  }
393 
394  // Store cosine of simulation cone radius and cone centre
395  m_mc_centre = dir;
397 
398  // Compute fraction of sky covered by the simulation cone
399  double fraction = 0.5 * (1.0 - m_mc_cos_radius);
400 
401  // Multiply by normalization value of model
402  double norm = fraction * value();
403 
404  // Return normalization
405  return (norm);
406 }
407 
408 
409 /***********************************************************************//**
410  * @brief Read model from XML element
411  *
412  * @param[in] xml XML element.
413  *
414  * Read the isotropic source model information from an XML element. The XML
415  * element is expected to have the following format:
416  *
417  * <spatialModel type="DiffuseIsotropic">
418  * <parameter name="Value" scale="1" value="1" min="1" max="1" free="0"/>
419  * </spatialModel>
420  *
421  ***************************************************************************/
423 {
424  // Get value parameter
425  const GXmlElement* par = gammalib::xml_get_par(G_READ, xml, m_value.name());
426 
427  // Read value parameter
428  m_value.read(*par);
429 
430  // Return
431  return;
432 }
433 
434 
435 /***********************************************************************//**
436  * @brief Write model into XML element
437  *
438  * @param[in] xml XML element.
439  *
440  * @exception GException::model_invalid_spatial
441  * Existing XML element is not of type "ConstantValue"
442  *
443  * Write the isotropic source model information into an XML element. The XML
444  * element will have the following format:
445  *
446  * <spatialModel type="DiffuseIsotropic">
447  * <parameter name="Value" scale="1" value="1" min="1" max="1" free="0"/>
448  * </spatialModel>
449  *
450  ***************************************************************************/
452 {
453  // Set model type
454  if (xml.attribute("type") == "") {
455  xml.attribute("type", type());
456  }
457 
458  // Verify model type
459  if (xml.attribute("type") != type()) {
461  "Spatial model is not of type \""+type()+"\".");
462  }
463 
464  // Get or create Value parameter
466 
467  // Write parameter
468  m_value.write(*par);
469 
470  // Return
471  return;
472 }
473 
474 
475 /***********************************************************************//**
476  * @brief Print isotropic source model information
477  *
478  * @param[in] chatter Chattiness (defaults to NORMAL).
479  * @return String containing model information.
480  ***************************************************************************/
481 std::string GModelSpatialDiffuseConst::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("=== GModelSpatialDiffuseConst ===");
491 
492  // Append parameters
493  result.append("\n"+gammalib::parformat("Number of parameters"));
494  result.append(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 model type
518  m_type = "DiffuseIsotropic";
519 
520  // Initialise Value
521  m_value.clear();
522  m_value.name("Value");
523  m_value.fix();
524  m_value.value(1.0);
525  m_value.scale(1.0);
526  m_value.range(0.0, 10.0);
527  m_value.gradient(0.0);
528  m_value.has_grad(true);
529 
530  // Initialise other members
531  m_mc_centre.clear();
532  m_mc_cos_radius = -1.0; // Set simulation cone to allsky
533 
534  // Set parameter pointer(s)
535  m_pars.clear();
536  m_pars.push_back(&m_value);
537 
538  // Initialise other members
539  m_region.clear();
540 
541  // Return
542  return;
543 }
544 
545 
546 /***********************************************************************//**
547  * @brief Copy class members
548  *
549  * @param[in] model Isotropic spatial model.
550  ***************************************************************************/
552 {
553  // Copy members
554  m_type = model.m_type;
555  m_value = model.m_value;
556  m_region = model.m_region;
557  m_mc_centre = model.m_mc_centre;
559 
560  // Set parameter pointer(s)
561  m_pars.clear();
562  m_pars.push_back(&m_value);
563 
564  // Return
565  return;
566 }
567 
568 
569 /***********************************************************************//**
570  * @brief Delete class members
571  ***************************************************************************/
573 {
574  // Return
575  return;
576 }
577 
578 
579 /***********************************************************************//**
580  * @brief Set boundary sky region
581  ***************************************************************************/
583 {
584  // Set sky region centre to (0,0)
585  m_region.centre(0.0, 0.0);
586 
587  // Set sky region radius to 180 degrees (all points included)
588  m_region.radius(180.0);
589 
590  // Return
591  return;
592 }
const double & factor_gradient(void) const
Return parameter gradient factor.
void init_members(void)
Initialise class members.
double norm(const GVector &vector)
Computes vector norm.
Definition: GVector.cpp:821
const std::string & name(void) const
Return parameter name.
GSkyRegionCircle m_region
Bounding circle.
double gradient(void) const
Return parameter gradient.
int size(void) const
Return number of parameters.
void free_members(void)
Delete 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
XML element node class.
Definition: GXmlElement.hpp:47
double acos(const double &arg)
Computes acos by avoiding NaN due to rounding errors.
Definition: GMath.cpp:69
Random number generator class.
Definition: GRan.hpp:44
Time class.
Definition: GTime.hpp:54
Gammalib tools definition.
Interface definition for the spatial model registry class.
bool is_free(void) const
Signal if parameter is free.
virtual std::string print(const GChatter &chatter=NORMAL) const
Print isotropic source model information.
void set_region(void) const
Set boundary sky region.
virtual void read(const GXmlElement &xml)
Read model from XML element.
const double & radius(void) const
Return circular region radius (in degrees)
virtual void write(GXmlElement &xml) const
Write model into XML element.
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.
void init_members(void)
Initialise class members.
const GXmlAttribute * attribute(const int &index) const
Return attribute.
virtual double mc_norm(const GSkyDir &dir, const double &radius) const
Return normalization of constant diffuse source for Monte Carlo simulations.
virtual ~GModelSpatialDiffuseConst(void)
Destructor.
std::vector< GModelPar * > m_pars
Parameter pointers.
void rotate_deg(const double &phi, const double &theta)
Rotate sky direction by zenith and azimuth angle.
Definition: GSkyDir.cpp:321
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
#define G_WRITE
double uniform(void)
Returns random double precision floating value in range 0 to 1.
Definition: GRan.cpp:242
#define G_READ
void clear(void)
Clear parameter.
#define G_MC_NORM
Spatial model registry class definition.
virtual void clear(void)
Clear isotropic spatial model.
GSkyDir m_mc_centre
Simulation cone centre.
GChatter
Definition: GTypemaps.hpp:33
double m_mc_cos_radius
Cosine of sim. cone radius.
Isotropic spatial model class interface definition.
virtual GModelSpatialDiffuse & operator=(const GModelSpatialDiffuse &model)
Assignment operator.
GModelSpatialDiffuseConst(void)
Void constructor.
void free_members(void)
Delete class members.
void copy_members(const GModelSpatialDiffuseConst &model)
Copy class members.
const GSkyDir & centre(void) const
Return circular region centre.
void clear(void)
Clear sky direction.
Definition: GSkyDir.cpp:142
void free_members(void)
Delete class members.
virtual std::string type(void) const
Return spatial model type.
virtual double eval(const GPhoton &photon, const bool &gradients=false) const
Evaluate isotropic spatial model value.
virtual GModelSpatialDiffuseConst * clone(void) const
Clone isotropic spatial model.
void range(const double &min, const double &max)
Set minimum and maximum parameter boundaries.
void read(const GXmlElement &xml)
Extract parameter attributes from XML element.
Definition: GModelPar.cpp:229
double value(void) const
Return parameter value.
const GModelSpatialDiffuseConst g_spatial_const_seed
GModelSpatialDiffuseConst & operator=(const GModelSpatialDiffuseConst &model)
Assignment operator.
Exception handler interface definition.
double value(void) const
Get model value.
const double fourpi
Definition: GMath.hpp:37
void init_members(void)
Initialise class members.
Abstract diffuse spatial model base class.
const double rad2deg
Definition: GMath.hpp:44
std::string parformat(const std::string &s, const int &indent=0)
Convert string in parameter format.
Definition: GTools.cpp:1022
Sky direction class.
Definition: GSkyDir.hpp:62
virtual GSkyDir mc(const GEnergy &energy, const GTime &time, GRan &ran) const
Return MC sky direction.
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
void autoscale(void)
Autoscale parameters.
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