GammaLib  1.7.0.dev
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
GModelSpectralConst.cpp
Go to the documentation of this file.
1 /***************************************************************************
2  * GModelSpectralConst.cpp - Spectral constant model class *
3  * ----------------------------------------------------------------------- *
4  * copyright (C) 2009-2016 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 GModelSpectralConst.cpp
23  * @brief Constant spectral 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 "GRan.hpp"
34 #include "GModelSpectralConst.hpp"
36 
37 /* __ Constants __________________________________________________________ */
38 
39 /* __ Globals ____________________________________________________________ */
40 const GModelSpectralConst g_spectral_const_seed1("Constant", "Normalization");
41 const GModelSpectralRegistry g_spectral_const_registry1(&g_spectral_const_seed1);
42 #if defined(G_LEGACY_XML_FORMAT)
43 const GModelSpectralConst g_spectral_const_seed2("ConstantValue", "Value");
44 const GModelSpectralRegistry g_spectral_const_registry2(&g_spectral_const_seed2);
45 #endif
46 
47 /* __ Method name definitions ____________________________________________ */
48 #define G_FLUX "GModelSpectralConst::flux(GEnergy&, GEnergy&)"
49 #define G_EFLUX "GModelSpectralConst::eflux(GEnergy&, GEnergy&)"
50 #define G_MC "GModelSpectralConst::mc(GEnergy&, GEnergy&, GTime&, GRan&)"
51 #define G_READ "GModelSpectralConst::read(GXmlElement&)"
52 #define G_WRITE "GModelSpectralConst::write(GXmlElement&)"
53 
54 /* __ Macros _____________________________________________________________ */
55 
56 /* __ Coding definitions _________________________________________________ */
57 
58 /* __ Debug definitions __________________________________________________ */
59 
60 
61 /*==========================================================================
62  = =
63  = Constructors/destructors =
64  = =
65  ==========================================================================*/
66 
67 /***********************************************************************//**
68  * @brief Void constructor
69  ***************************************************************************/
71 {
72  // Initialise members
73  init_members();
74 
75  // Return
76  return;
77 }
78 
79 
80 /***********************************************************************//**
81  * @brief Model type and parameter name constructor
82  *
83  * @param[in] type Model type.
84  * @param[in] value Name of value parameter.
85  ***************************************************************************/
87  const std::string& value) :
89 {
90  // Initialise members
91  init_members();
92 
93  // Set model type
94  m_type = type;
95 
96  // Set parameter names
97  m_norm.name(value);
98 
99  // Return
100  return;
101 }
102 
103 
104 /***********************************************************************//**
105  * @brief XML constructor
106  *
107  * @param[in] xml XML element.
108  *
109  * Constructs constant spectral 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 Model value (ph/cm2/s/MeV).
131  *
132  * Constructs constant spectral model by setting the model value.
133  ***************************************************************************/
136 {
137  // Initialise members
138  init_members();
139 
140  // Set model value
141  m_norm.value(value);
142 
143  // Return
144  return;
145 }
146 
147 
148 /***********************************************************************//**
149  * @brief Copy constructor
150  *
151  * @param[in] model Spectral constant model.
152  ***************************************************************************/
154  GModelSpectral(model)
155 {
156  // Initialise members
157  init_members();
158 
159  // Copy members
160  copy_members(model);
161 
162  // Return
163  return;
164 }
165 
166 
167 /***********************************************************************//**
168  * @brief Destructor
169  ***************************************************************************/
171 {
172  // Free members
173  free_members();
174 
175  // Return
176  return;
177 }
178 
179 
180 /*==========================================================================
181  = =
182  = Operators =
183  = =
184  ==========================================================================*/
185 
186 /***********************************************************************//**
187  * @brief Assignment operator
188  *
189  * @param[in] model Constant spectral model.
190  * @return Constant spectral model.
191  ***************************************************************************/
193 {
194  // Execute only if object is not identical
195  if (this != &model) {
196 
197  // Copy base class members
198  this->GModelSpectral::operator=(model);
199 
200  // Free members
201  free_members();
202 
203  // Initialise members
204  init_members();
205 
206  // Copy members
207  copy_members(model);
208 
209  } // endif: object was not identical
210 
211  // Return
212  return *this;
213 }
214 
215 
216 /*==========================================================================
217  = =
218  = Public methods =
219  = =
220  ==========================================================================*/
221 
222 /***********************************************************************//**
223  * @brief Clear constant spectral model
224  ***************************************************************************/
226 {
227  // Free class members (base and derived classes, derived class first)
228  free_members();
230 
231  // Initialise members
233  init_members();
234 
235  // Return
236  return;
237 }
238 
239 
240 /***********************************************************************//**
241  * @brief Clone constant spectral model
242  *
243  * @return Pointer to deep copy of constant spectral model.
244  ***************************************************************************/
246 {
247  // Clone constant spectral model
248  return new GModelSpectralConst(*this);
249 }
250 
251 
252 /***********************************************************************//**
253  * @brief Evaluate model value
254  *
255  * @param[in] srcEng True photon energy.
256  * @param[in] srcTime True photon arrival time.
257  * @param[in] gradients Compute gradients?
258  * @return Model value (ph/cm2/s/MeV).
259  *
260  * Evaluates
261  *
262  * \f[
263  * S_{\rm E}(E | t) = {\tt m\_norm}
264  * \f]
265  *
266  * where
267  * \f${\tt m\_norm}\f$ is the normalization constant in units of
268  * ph/cm2/s/MeV.
269  *
270  * If the @p gradients flag is true the method will also compute the
271  * partial derivatives of the model with respect to the parameters using
272  *
273  * \f[
274  * \frac{\delta S_{\rm E}(E | t)}{\delta {\tt m\_norm}} = 1
275  * \f]
276  ***************************************************************************/
277 double GModelSpectralConst::eval(const GEnergy& srcEng,
278  const GTime& srcTime,
279  const bool& gradients) const
280 {
281  // Compute function value
282  double value = m_norm.value();
283 
284  // Optionally compute gradients
285  if (gradients) {
286 
287  // Compute partial derivatives of the parameter values
288  double g_norm = (m_norm.is_free()) ? m_norm.scale() : 0.0;
289 
290  // Set factor gradient
291  m_norm.factor_gradient(g_norm);
292 
293  } // endif: gradient computation was requested
294 
295  // Return
296  return value;
297 }
298 
299 
300 /***********************************************************************//**
301  * @brief Returns model photon flux between [emin, emax] (ph/cm2/s)
302  *
303  * @param[in] emin Minimum photon energy.
304  * @param[in] emax Maximum photon energy.
305  * @return Photon flux (ph/cm2/s).
306  *
307  * Computes
308  *
309  * \f[
310  * \int_{\tt emin}^{\tt emax} S_{\rm E}(E | t) dE
311  * \f]
312  *
313  * where
314  * - [@p emin, @p emax] is an energy interval, and
315  * - \f$S_{\rm E}(E | t)\f$ is the spectral model (ph/cm2/s/MeV).
316  * The integration is done analytically.
317  ***************************************************************************/
319  const GEnergy& emax) const
320 {
321  // Initialise flux
322  double flux = 0.0;
323 
324  // Compute only if integration range is valid
325  if (emin < emax) {
326 
327  // Compute flux for a constant model
328  flux = m_norm.value() * (emax.MeV() - emin.MeV());
329 
330  } // endif: integration range was valid
331 
332  // Return
333  return flux;
334 }
335 
336 
337 /***********************************************************************//**
338  * @brief Returns model energy flux between [emin, emax] (erg/cm2/s)
339  *
340  * @param[in] emin Minimum photon energy.
341  * @param[in] emax Maximum photon energy.
342  * @return Energy flux (erg/cm2/s).
343  *
344  * Computes
345  *
346  * \f[
347  * \int_{\tt emin}^{\tt emax} S_{\rm E}(E | t) E \, dE
348  * \f]
349  *
350  * where
351  * - [@p emin, @p emax] is an energy interval, and
352  * - \f$S_{\rm E}(E | t)\f$ is the spectral model (ph/cm2/s/MeV).
353  * The integration is done analytically.
354  ***************************************************************************/
356  const GEnergy& emax) const
357 {
358  // Initialise flux
359  double eflux = 0.0;
360 
361  // Compute only if integration range is valid
362  if (emin < emax) {
363 
364  // Compute flux for a constant model
365  eflux = m_norm.value() * 0.5 * (emax.MeV()*emax.MeV() -
366  emin.MeV()*emin.MeV());
367 
368  // Convert from MeV/cm2/s to erg/cm2/s
369  eflux *= gammalib::MeV2erg;
370 
371  } // endif: integration range was valid
372 
373  // Return
374  return eflux;
375 }
376 
377 
378 /***********************************************************************//**
379  * @brief Returns MC energy between [emin, emax]
380  *
381  * @param[in] emin Minimum photon energy.
382  * @param[in] emax Maximum photon energy.
383  * @param[in] time True photon arrival time.
384  * @param[in,out] ran Random number generator.
385  * @return Energy.
386  *
387  * @exception GException::erange_invalid
388  * Energy range is invalid (emin < emax required).
389  *
390  * Returns Monte Carlo energy by randomly drawing from a constant between
391  * the minimum and maximum photon energy.
392  ***************************************************************************/
394  const GEnergy& emax,
395  const GTime& time,
396  GRan& ran) const
397 {
398  // Throw an exception if energy range is invalid
399  if (emin >= emax) {
400  throw GException::erange_invalid(G_MC, emin.MeV(), emax.MeV(),
401  "Minimum energy < maximum energy required.");
402  }
403 
404  // Allocate energy
405  GEnergy energy;
406 
407  // Get uniform value between 0 and 1
408  double u = ran.uniform();
409 
410  // Map into [emin,emax] range
411  energy.MeV((emax.MeV() - emin.MeV()) * u + emin.MeV());
412 
413  // Return energy
414  return energy;
415 }
416 
417 
418 /***********************************************************************//**
419  * @brief Read model from XML element
420  *
421  * @param[in] xml XML element.
422  *
423  * Reads the spectral information from an XML element.
424  ***************************************************************************/
426 {
427  // Get parameter pointers
429 
430  // Read parameters
431  m_norm.read(*norm);
432 
433  // Return
434  return;
435 }
436 
437 
438 /***********************************************************************//**
439  * @brief Write model into XML element
440  *
441  * @param[in] xml XML element.
442  *
443  * @exception GException::model_invalid_spectral
444  * Existing XML element is not of type "ConstantValue"
445  *
446  * Writes the spectral information into an XML element.
447  ***************************************************************************/
449 {
450  // Set model type
451  if (xml.attribute("type") == "") {
452  xml.attribute("type", type());
453  }
454 
455  // Verify model type
456  if (xml.attribute("type") != type()) {
458  "Spectral model is not of type \""+type()+"\".");
459  }
460 
461  // Get normalisation parameter
463 
464  // Write normalisation parameter
465  m_norm.write(*par);
466 
467  // Return
468  return;
469 }
470 
471 
472 /***********************************************************************//**
473  * @brief Print spectral model information
474  *
475  * @param[in] chatter Chattiness (defaults to NORMAL).
476  * @return String containing spectral model information.
477  ***************************************************************************/
478 std::string GModelSpectralConst::print(const GChatter& chatter) const
479 {
480  // Initialise result string
481  std::string result;
482 
483  // Continue only if chatter is not silent
484  if (chatter != SILENT) {
485 
486  // Append header
487  result.append("=== GModelSpectralConst ===");
488 
489  // Append model content
490  result.append("\n"+gammalib::parformat("Number of parameters"));
491  result.append(gammalib::str(size()));
492  for (int i = 0; i < size(); ++i) {
493  result.append("\n"+m_pars[i]->print(chatter));
494  }
495 
496  } // endif: chatter was not silent
497 
498  // Return result
499  return result;
500 }
501 
502 
503 /*==========================================================================
504  = =
505  = Private methods =
506  = =
507  ==========================================================================*/
508 
509 /***********************************************************************//**
510  * @brief Initialise class members
511  ***************************************************************************/
513 {
514  // Initialise model type
515  m_type = "Constant";
516 
517  // Initialise constant normalisation
518  m_norm.clear();
519  m_norm.name("Normalization");
520  m_norm.scale(1.0);
521  m_norm.value(1.0); // default: 1.0
522  m_norm.range(0.0, 1000.0); // range: [0, 1000]
523  m_norm.free();
524  m_norm.gradient(0.0);
525  m_norm.has_grad(true);
526 
527  // Set parameter pointer(s)
528  m_pars.clear();
529  m_pars.push_back(&m_norm);
530 
531  // Return
532  return;
533 }
534 
535 
536 /***********************************************************************//**
537  * @brief Copy class members
538  *
539  * @param[in] model Spectral constant model.
540  ***************************************************************************/
542 {
543  // Copy members
544  m_type = model.m_type;
545  m_norm = model.m_norm;
546 
547  // Set parameter pointer(s)
548  m_pars.clear();
549  m_pars.push_back(&m_norm);
550 
551  // Return
552  return;
553 }
554 
555 
556 /***********************************************************************//**
557  * @brief Delete class members
558  ***************************************************************************/
560 {
561  // Return
562  return;
563 }
virtual void read(const GXmlElement &xml)
Read model from XML element.
void init_members(void)
Initialise class members.
const double & factor_gradient(void) const
Return parameter gradient factor.
double norm(const GVector &vector)
Computes vector norm.
Definition: GVector.cpp:821
virtual GModelSpectralConst & operator=(const GModelSpectralConst &model)
Assignment operator.
const std::string & name(void) const
Return parameter name.
Random number generator class definition.
Abstract spectral model base class.
virtual GModelSpectralConst * clone(void) const
Clone constant spectral model.
double gradient(void) const
Return parameter gradient.
virtual double eval(const GEnergy &srcEng, const GTime &srcTime=GTime(), const bool &gradients=false) const
Evaluate model value.
int size(void) const
Return number of parameters.
#define G_MC
virtual GModelSpectral & operator=(const GModelSpectral &model)
Assignment operator.
void write(GXmlElement &xml) const
Set or update parameter attributes in XML element.
Definition: GModelPar.cpp:353
std::vector< GModelPar * > m_pars
Parameter pointers.
virtual double eflux(const GEnergy &emin, const GEnergy &emax) const
Returns model energy flux between emin, emax
std::string m_type
Model type.
XML element node class.
Definition: GXmlElement.hpp:47
double value(void) const
Return model value.
Spectral model registry class definition.
const GModelSpectralConst g_spectral_const_seed1("Constant","Normalization")
Random number generator class.
Definition: GRan.hpp:44
double MeV(void) const
Return energy in MeV.
Definition: GEnergy.cpp:321
Time class.
Definition: GTime.hpp:54
Gammalib tools definition.
void copy_members(const GModelSpectralConst &model)
Copy class members.
bool is_free(void) const
Signal if parameter is free.
const double & scale(void) const
Return parameter scale.
GModelSpectralConst(void)
Void constructor.
virtual void clear(void)
Clear constant spectral model.
bool has_grad(void) const
Signal if parameter gradient is computed analytically.
const GXmlAttribute * attribute(const int &index) const
Return attribute.
void free(void)
Free 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
Constant spectral model class interface definition.
double uniform(void)
Returns random double precision floating value in range 0 to 1.
Definition: GRan.cpp:242
void clear(void)
Clear parameter.
#define G_READ
GChatter
Definition: GTypemaps.hpp:33
virtual double flux(const GEnergy &emin, const GEnergy &emax) const
Returns model photon flux between emin, emax
virtual std::string print(const GChatter &chatter=NORMAL) const
Print spectral model information.
Interface definition for the spectral model registry class.
virtual std::string type(void) const
Return model type.
void init_members(void)
Initialise class members.
const double MeV2erg
Definition: GTools.hpp:45
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.
#define G_WRITE
Exception handler interface definition.
virtual ~GModelSpectralConst(void)
Destructor.
GModelPar m_norm
Normalization factor.
virtual void write(GXmlElement &xml) const
Write model into XML element.
std::string parformat(const std::string &s, const int &indent=0)
Convert string in parameter format.
Definition: GTools.cpp:1022
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 free_members(void)
Delete class members.
Constant spectral model class.
void free_members(void)
Delete class members.
virtual GEnergy mc(const GEnergy &emin, const GEnergy &emax, const GTime &time, GRan &ran) const
Returns MC energy between [emin, emax].
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