GammaLib  1.7.0.dev
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
GModelSpectralFunc.cpp
Go to the documentation of this file.
1 /***************************************************************************
2  * GModelSpectralFunc.cpp - Spectral function model class *
3  * ----------------------------------------------------------------------- *
4  * copyright (C) 2010-2017 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 GModelSpectralFunc.cpp
23  * @brief Spectral function 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 "GCsv.hpp"
35 #include "GRan.hpp"
36 #include "GModelSpectralFunc.hpp"
38 
39 /* __ Constants __________________________________________________________ */
40 
41 /* __ Globals ____________________________________________________________ */
43 const GModelSpectralRegistry g_spectral_func_registry(&g_spectral_func_seed);
44 
45 /* __ Method name definitions ____________________________________________ */
46 #define G_FLUX "GModelSpectralFunc::flux(GEnergy&, GEnergy&)"
47 #define G_EFLUX "GModelSpectralFunc::eflux(GEnergy&, GEnergy&)"
48 #define G_MC "GModelSpectralFunc::mc(GEnergy&, GEnergy&, GTime&, GRan&)"
49 #define G_READ "GModelSpectralFunc::read(GXmlElement&)"
50 #define G_WRITE "GModelSpectralFunc::write(GXmlElement&)"
51 #define G_LOAD_NODES "GModelSpectralFunc::load_nodes(GFilename&)"
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  ***************************************************************************/
70 {
71  // Initialise members
72  init_members();
73 
74  // Return
75  return;
76 }
77 
78 
79 /***********************************************************************//**
80  * @brief File constructor
81  *
82  * @param[in] filename File name of nodes.
83  * @param[in] norm Normalization factor.
84  *
85  * Constructs spectral file function model from a list of nodes that is found
86  * in the specified ASCII file. See the load_nodes() method for more
87  * information about the expected structure of the file.
88  ***************************************************************************/
90  const double& norm) :
92 {
93  // Initialise members
94  init_members();
95 
96  // Load nodes
97  load_nodes(filename);
98 
99  // Set normalization
100  m_norm.value(norm);
101 
102  // Return
103  return;
104 }
105 
106 
107 /***********************************************************************//**
108  * @brief XML constructor
109  *
110  * @param[in] xml XML element.
111  *
112  * Constructs spectral file function model by extracting information from an
113  * XML element. See the read() method for more information about the expected
114  * structure of the XML element.
115  ***************************************************************************/
118 {
119  // Initialise members
120  init_members();
121 
122  // Read information from XML element
123  read(xml);
124 
125  // Return
126  return;
127 }
128 
129 
130 /***********************************************************************//**
131  * @brief Copy constructor
132  *
133  * @param[in] model File function model.
134  ***************************************************************************/
136  GModelSpectral(model)
137 {
138  // Initialise members
139  init_members();
140 
141  // Copy members
142  copy_members(model);
143 
144  // Return
145  return;
146 }
147 
148 
149 /***********************************************************************//**
150  * @brief Destructor
151  ***************************************************************************/
153 {
154  // Free members
155  free_members();
156 
157  // Return
158  return;
159 }
160 
161 
162 /*==========================================================================
163  = =
164  = Operators =
165  = =
166  ==========================================================================*/
167 
168 /***********************************************************************//**
169  * @brief Assignment operator
170  *
171  * @param[in] model File function model.
172  * @return File function model.
173  ***************************************************************************/
175 {
176  // Execute only if object is not identical
177  if (this != &model) {
178 
179  // Copy base class members
180  this->GModelSpectral::operator=(model);
181 
182  // Free members
183  free_members();
184 
185  // Initialise members
186  init_members();
187 
188  // Copy members
189  copy_members(model);
190 
191  } // endif: object was not identical
192 
193  // Return
194  return *this;
195 }
196 
197 
198 /*==========================================================================
199  = =
200  = Public methods =
201  = =
202  ==========================================================================*/
203 
204 /***********************************************************************//**
205  * @brief Clear file function
206 ***************************************************************************/
208 {
209  // Free class members (base and derived classes, derived class first)
210  free_members();
212 
213  // Initialise members
215  init_members();
216 
217  // Return
218  return;
219 }
220 
221 
222 /***********************************************************************//**
223  * @brief Clone file function
224 ***************************************************************************/
226 {
227  // Clone file function
228  return new GModelSpectralFunc(*this);
229 }
230 
231 
232 /***********************************************************************//**
233  * @brief Evaluate function
234  *
235  * @param[in] srcEng True photon energy.
236  * @param[in] srcTime True photon arrival time.
237  * @param[in] gradients Compute gradients?
238  * @return Model value (ph/cm2/s/MeV).
239  *
240  * Evaluates
241  *
242  * \f[
243  * S_{\rm E}(E | t) = {\tt m\_norm}
244  * \f]
245  *
246  * where
247  * - \f${\tt m\_norm}\f$ is the normalization factor.
248  *
249  * If the @p gradients flag is true the method will also compute the
250  * partial derivatives of the model with respect to the parameters using
251  *
252  * \f[
253  * \frac{\delta S_{\rm E}(E | t)}{\delta {\tt m\_norm}} =
254  * \frac{S_{\rm E}(E | t)}{{\tt m\_norm}}
255  * \f]
256  ***************************************************************************/
257 double GModelSpectralFunc::eval(const GEnergy& srcEng,
258  const GTime& srcTime,
259  const bool& gradients) const
260 {
261  // Interpolate function. This is done in log10-log10 space, but the
262  // linear value is returned.
263  double arg = m_log_nodes.interpolate(srcEng.log10MeV(), m_log_values);
264  double func = std::pow(10.0, arg);
265 
266  // Compute function value
267  double value = m_norm.value() * func;
268 
269  // Optionally compute gradients
270  if (gradients) {
271 
272  // Compute partial derivatives of the parameter values
273  double g_norm = (m_norm.is_free()) ? m_norm.scale() * func : 0.0;
274 
275  // Set gradients
276  m_norm.factor_gradient(g_norm);
277 
278  } // endif: gradient computation was requested
279 
280  // Compile option: Check for NaN/Inf
281  #if defined(G_NAN_CHECK)
282  if (gammalib::is_notanumber(value) || gammalib::is_infinite(value)) {
283  std::cout << "*** ERROR: GModelSpectralFunc::eval";
284  std::cout << "(srcEng=" << srcEng;
285  std::cout << ", srcTime=" << srcTime << "):";
286  std::cout << " NaN/Inf encountered";
287  std::cout << " (value=" << value;
288  std::cout << ", norm=" << norm();
289  std::cout << ", func=" << func;
290  std::cout << ")" << std::endl;
291  }
292  #endif
293 
294  // Return
295  return value;
296 }
297 
298 
299 /***********************************************************************//**
300  * @brief Returns model photon flux between [emin, emax] (units: ph/cm2/s)
301  *
302  * @param[in] emin Minimum photon energy.
303  * @param[in] emax Maximum photon energy.
304  * @return Photon flux (ph/cm2/s).
305  *
306  * Computes
307  *
308  * \f[
309  * \int_{\tt emin}^{\tt emax} S_{\rm E}(E | t) dE
310  * \f]
311  *
312  * where
313  * - [@p emin, @p emax] is an energy interval, and
314  * - \f$S_{\rm E}(E | t)\f$ is the spectral model (ph/cm2/s/MeV).
315  ***************************************************************************/
316 double GModelSpectralFunc::flux(const GEnergy& emin, const GEnergy& emax) const
317 {
318  // Initialise flux
319  double flux = 0.0;
320 
321  // Compute only if integration range is valid
322  if (emin < emax) {
323 
324  // Get energy range in MeV
325  double e_min = emin.MeV();
326  double e_max = emax.MeV();
327 
328  // Determine left node index for minimum energy
329  m_lin_nodes.set_value(e_min);
330  int inx_emin = m_lin_nodes.inx_left();
331 
332  // Determine left node index for maximum energy
333  m_lin_nodes.set_value(e_max);
334  int inx_emax = m_lin_nodes.inx_left();
335 
336  // If both energies are within the same nodes then simply
337  // integrate over the energy interval using the appropriate power
338  // law parameters
339  if (inx_emin == inx_emax) {
340  flux = m_prefactor[inx_emin] *
342  e_max,
343  m_epivot[inx_emin],
344  m_gamma[inx_emin]);
345  }
346 
347  // ... otherwise integrate over the nodes where emin and emax
348  // resides and all the remaining nodes
349  else {
350 
351  // If we are requested to extrapolate beyond first node,
352  // the use the first nodes lower energy and upper energy
353  // boundary for the initial integration.
354  int i_start = (e_min < m_lin_nodes[0]) ? inx_emin : inx_emin+1;
355 
356  // Integrate from emin to the node boundary
357  flux = m_prefactor[inx_emin] *
359  m_lin_nodes[i_start],
360  m_epivot[inx_emin],
361  m_gamma[inx_emin]);
362 
363  // Integrate over all nodes between
364  for (int i = i_start; i < inx_emax; ++i) {
365  flux += m_flux[i];
366  }
367 
368  // Integrate from node boundary to emax
369  flux += m_prefactor[inx_emax] *
371  e_max,
372  m_epivot[inx_emax],
373  m_gamma[inx_emax]);
374 
375  } // endelse: emin and emax not between same nodes
376 
377  // Multiply flux by normalisation factor
378  flux *= m_norm.value();
379 
380  } // endif: integration range was valid
381 
382  // Return
383  return flux;
384 }
385 
386 
387 /***********************************************************************//**
388  * @brief Returns model energy flux between [emin, emax] (units: erg/cm2/s)
389  *
390  * @param[in] emin Minimum photon energy.
391  * @param[in] emax Maximum photon energy.
392  * @return Energy flux (erg/cm2/s).
393  *
394  * Computes
395  *
396  * \f[
397  * \int_{\tt emin}^{\tt emax} S_{\rm E}(E | t) E \, dE
398  * \f]
399  *
400  * where
401  * - [@p emin, @p emax] is an energy interval, and
402  * - \f$S_{\rm E}(E | t)\f$ is the spectral model (ph/cm2/s/MeV).
403  ***************************************************************************/
404 double GModelSpectralFunc::eflux(const GEnergy& emin, const GEnergy& emax) const
405 {
406  // Initialise flux
407  double eflux = 0.0;
408 
409  // Compute only if integration range is valid
410  if (emin < emax) {
411 
412  // Get energy range in MeV
413  double e_min = emin.MeV();
414  double e_max = emax.MeV();
415 
416  // Determine left node index for minimum energy
417  m_lin_nodes.set_value(e_min);
418  int inx_emin = m_lin_nodes.inx_left();
419 
420  // Determine left node index for maximum energy
421  m_lin_nodes.set_value(e_max);
422  int inx_emax = m_lin_nodes.inx_left();
423 
424  // If both energies are within the same nodes then simply
425  // integrate over the energy interval using the appropriate power
426  // law parameters
427  if (inx_emin == inx_emax) {
428  eflux = m_prefactor[inx_emin] *
430  e_max,
431  m_epivot[inx_emin],
432  m_gamma[inx_emin]) *
434  }
435 
436  // ... otherwise integrate over the nodes where emin and emax
437  // resides and all the remaining nodes
438  else {
439 
440  // If we are requested to extrapolate beyond first node,
441  // the use the first nodes lower energy and upper energy
442  // boundary for the initial integration.
443  int i_start = (e_min < m_lin_nodes[0]) ? inx_emin : inx_emin+1;
444 
445  // Integrate from emin to the node boundary
446  eflux = m_prefactor[inx_emin] *
448  m_lin_nodes[i_start],
449  m_epivot[inx_emin],
450  m_gamma[inx_emin]) *
452 
453  // Integrate over all nodes between
454  for (int i = i_start; i < inx_emax; ++i) {
455  eflux += m_eflux[i];
456  }
457 
458  // Integrate from node boundary to emax
459  eflux += m_prefactor[inx_emax] *
461  e_max,
462  m_epivot[inx_emax],
463  m_gamma[inx_emax]) *
465 
466  } // endelse: emin and emax not between same nodes
467 
468  // Multiply flux by normalisation factor
469  eflux *= norm();
470 
471  } // endif: integration range was valid
472 
473  // Return flux
474  return eflux;
475 }
476 
477 
478 /***********************************************************************//**
479  * @brief Returns MC energy between [emin, emax]
480  *
481  * @param[in] emin Minimum photon energy.
482  * @param[in] emax Maximum photon energy.
483  * @param[in] time True photon arrival time.
484  * @param[in,out] ran Random number generator.
485  * @return Energy.
486  *
487  * @exception GException::erange_invalid
488  * Energy range is invalid (emin < emax required).
489  *
490  * Returns Monte Carlo energy by randomly drawing from a broken power law
491  * defined by the file function.
492  ***************************************************************************/
494  const GEnergy& emax,
495  const GTime& time,
496  GRan& ran) const
497 {
498  // Throw an exception if energy range is invalid
499  if (emin >= emax) {
500  throw GException::erange_invalid(G_MC, emin.MeV(), emax.MeV(),
501  "Minimum energy < maximum energy required.");
502  }
503 
504  // Allocate energy
505  GEnergy energy;
506 
507  // Update cache
508  mc_update(emin, emax);
509 
510  // Determine in which bin we reside
511  int inx = 0;
512  if (m_mc_cum.size() > 1) {
513  double u = ran.uniform();
514  for (inx = m_mc_cum.size()-1; inx > 0; --inx) {
515  if (m_mc_cum[inx-1] <= u) {
516  break;
517  }
518  }
519  }
520 
521  // Get random energy for specific bin
522  if (m_mc_exp[inx] != 0.0) {
523  double e_min = m_mc_min[inx];
524  double e_max = m_mc_max[inx];
525  double u = ran.uniform();
526  double eng = (u > 0.0)
527  ? std::exp(std::log(u * (e_max - e_min) + e_min) / m_mc_exp[inx])
528  : 0.0;
529  energy.MeV(eng);
530  }
531  else {
532  double e_min = m_mc_min[inx];
533  double e_max = m_mc_max[inx];
534  double u = ran.uniform();
535  double eng = std::exp(u * (e_max - e_min) + e_min);
536  energy.MeV(eng);
537  }
538 
539  // Return energy
540  return energy;
541 }
542 
543 
544 /***********************************************************************//**
545  * @brief Read model from XML element
546  *
547  * @param[in] xml XML element containing power law model information.
548  *
549  * @exception GException::model_invalid_parnum
550  * Invalid number of model parameters found in XML element.
551  * @exception GException::model_invalid_parnames
552  * Invalid model parameter name found in XML element.
553  *
554  * Reads the spectral information from an XML element. The format of the XML
555  * elements is
556  *
557  * <spectrum type="FileFunction" file="..">
558  * <parameter name="Normalization" scale=".." value=".." min=".." max=".." free=".."/>
559  * </spectrum>
560  ***************************************************************************/
562 {
563  // Verify that XML element has exactly 1 parameter
564  if (xml.elements() != 1 || xml.elements("parameter") != 1) {
566  "Spectral function requires exactly 1 parameter.");
567  }
568 
569  // Get parameter element
570  const GXmlElement* par = xml.element("parameter", 0);
571 
572  // Get value
573  if (par->attribute("name") == "Normalization") {
574  m_norm.read(*par);
575  }
576  else {
578  "Require \"Normalization\" parameter.");
579  }
580 
581  // Load nodes from file
583 
584  // Return
585  return;
586 }
587 
588 
589 /***********************************************************************//**
590  * @brief Write model into XML element
591  *
592  * @param[in] xml XML element into which model information is written.
593  *
594  * @exception GException::model_invalid_spectral
595  * Existing XML element is not of type "FileFunction"
596  * @exception GException::model_invalid_parnum
597  * Invalid number of model parameters or nodes found in XML element.
598  * @exception GException::model_invalid_parnames
599  * Invalid model parameter names found in XML element.
600  *
601  * Writes the spectral information into an XML element. The format of the XML
602  * element is
603  *
604  * <spectrum type="FileFunction" file="..">
605  * <parameter name="Normalization" scale=".." value=".." min=".." max=".." free=".."/>
606  * </spectrum>
607  *
608  * Note that the function nodes will not be written since they will not be
609  * altered by any method.
610  ***************************************************************************/
612 {
613  // Set model type
614  if (xml.attribute("type") == "") {
615  xml.attribute("type", "FileFunction");
616  }
617 
618  // Verify model type
619  if (xml.attribute("type") != "FileFunction") {
621  "Spectral model is not of type \"FileFunction\".");
622  }
623 
624  // If XML element has 0 nodes then append 1 parameter node
625  if (xml.elements() == 0) {
626  xml.append(GXmlElement("parameter name=\"Normalization\""));
627  }
628 
629  // Verify that XML element has exactly 1 parameter
630  if (xml.elements() != 1 || xml.elements("parameter") != 1) {
632  "Spectral function requires exactly 1 parameter.");
633  }
634 
635  // Get parameter element
636  GXmlElement* par = xml.element("parameter", 0);
637 
638  // Set parameter
639  if (par->attribute("name") == "Normalization") {
640  m_norm.write(*par);
641  }
642  else {
644  "Require \"Normalization\" parameter.");
645  }
646 
647  // Set file attribute
649 
650  // Return
651  return;
652 }
653 
654 
655 /***********************************************************************//**
656  * @brief Print file function information
657  *
658  * @param[in] chatter Chattiness (defaults to NORMAL).
659  * @return String containing file function information.
660  ***************************************************************************/
661 std::string GModelSpectralFunc::print(const GChatter& chatter) const
662 {
663  // Initialise result string
664  std::string result;
665 
666  // Continue only if chatter is not silent
667  if (chatter != SILENT) {
668 
669  // Append header
670  result.append("=== GModelSpectralFunc ===");
671 
672  // Append information
673  result.append("\n"+gammalib::parformat("Function file"));
674  result.append(m_filename.url());
675  result.append("\n"+gammalib::parformat("Number of nodes"));
676  result.append(gammalib::str(m_lin_nodes.size()));
677  result.append("\n"+gammalib::parformat("Number of parameters"));
678  result.append(gammalib::str(size()));
679  for (int i = 0; i < size(); ++i) {
680  result.append("\n"+m_pars[i]->print(chatter));
681  }
682 
683  // Append node information
684  for (int i = 0; i < m_prefactor.size(); ++i) {
685  result.append("\n"+gammalib::parformat("Node "+gammalib::str(i+1)));
686  result.append("Epivot="+gammalib::str(m_epivot[i]));
687  result.append(" Prefactor="+gammalib::str(m_prefactor[i]));
688  result.append(" Gamma="+gammalib::str(m_gamma[i]));
689  result.append(" Flux="+gammalib::str(m_flux[i]));
690  result.append(" EFlux="+gammalib::str(m_eflux[i]));
691  }
692 
693  } // endif: chatter was not silent
694 
695  // Return result
696  return result;
697 }
698 
699 
700 /*==========================================================================
701  = =
702  = Private methods =
703  = =
704  ==========================================================================*/
705 
706 /***********************************************************************//**
707  * @brief Initialise class members
708  ***************************************************************************/
710 {
711  // Initialise powerlaw normalisation
712  m_norm.clear();
713  m_norm.name("Normalization");
714  m_norm.scale(1.0);
715  m_norm.value(1.0);
716  m_norm.range(0.0,1000.0);
717  m_norm.free();
718  m_norm.gradient(0.0);
719  m_norm.has_grad(true);
720 
721  // Set parameter pointer(s)
722  m_pars.clear();
723  m_pars.push_back(&m_norm);
724 
725  // Initialise members
726  m_lin_nodes.clear();
727  m_log_nodes.clear();
728  m_lin_values.clear();
729  m_log_values.clear();
730  m_filename.clear();
731 
732  // Initialise cache
733  m_mc_emin.clear();
734  m_mc_emax.clear();
735  m_mc_cum.clear();
736  m_mc_min.clear();
737  m_mc_max.clear();
738  m_mc_exp.clear();
739 
740  // Return
741  return;
742 }
743 
744 
745 /***********************************************************************//**
746  * @brief Copy class members
747  *
748  * @param[in] model Spectral function model.
749  ***************************************************************************/
751 {
752  // Copy members
753  m_norm = model.m_norm;
754  m_lin_nodes = model.m_lin_nodes;
755  m_log_nodes = model.m_log_nodes;
756  m_lin_values = model.m_lin_values;
757  m_log_values = model.m_log_values;
758  m_filename = model.m_filename;
759 
760  // Copy pre-computation cache
761  m_prefactor = model.m_prefactor;
762  m_gamma = model.m_gamma;
763  m_epivot = model.m_epivot;
764  m_flux = model.m_flux;
765  m_eflux = model.m_eflux;
766 
767  // Copy MC cache
768  m_mc_emin = model.m_mc_emin;
769  m_mc_emax = model.m_mc_emax;
770  m_mc_cum = model.m_mc_cum;
771  m_mc_min = model.m_mc_min;
772  m_mc_max = model.m_mc_max;
773  m_mc_exp = model.m_mc_exp;
774 
775  // Set parameter pointer(s)
776  m_pars.clear();
777  m_pars.push_back(&m_norm);
778 
779  // Return
780  return;
781 }
782 
783 
784 /***********************************************************************//**
785  * @brief Delete class members
786  ***************************************************************************/
788 {
789  // Return
790  return;
791 }
792 
793 
794 /***********************************************************************//**
795  * @brief Load nodes from file
796  *
797  * @param[in] filename File name.
798  *
799  * @exception GException::file_function_data
800  * File contains less than 2 nodes
801  * @exception GException::file_function_columns
802  * File contains less than 2 columns
803  * @exception GException::file_function_value
804  * File contains invalid value
805  *
806  * The file function is stored as a column separated value table (CSV) in an
807  * ASCII file with (at least) 2 columns. The first column specifies the
808  * energy in MeV while the second column specifies the intensity at this
809  * energy in units of ph/cm2/s/MeV.
810  * The node energies and values will be stored both linearly and as log10.
811  * The log10 storing requires that node energies and node values are
812  * positive. Also, at least 2 nodes and 2 columns are required in the file
813  * function.
814  ***************************************************************************/
816 {
817  // Clear nodes and values
818  m_lin_nodes.clear();
819  m_log_nodes.clear();
820  m_lin_values.clear();
821  m_log_values.clear();
822 
823  // Set filename
825 
826  // Continue only if filename is not empty
827  if (!filename.is_empty()) {
828 
829  // Load file
830  GCsv csv = GCsv(filename.url());
831 
832  // Check if there are at least 2 nodes
833  if (csv.nrows() < 2) {
835  csv.nrows());
836  }
837 
838  // Check if there are at least 2 columns
839  if (csv.ncols() < 2) {
841  csv.ncols());
842  }
843 
844  // Setup nodes
845  double last_energy = 0.0;
846  for (int i = 0; i < csv.nrows(); ++i) {
847 
848  // Get log10 of node energy and value. Make sure they are valid.
849  double log10energy;
850  double log10value;
851  if (csv.real(i,0) > 0) {
852  log10energy = std::log10(csv.real(i,0));
853  }
854  else {
856  csv.real(i,0), "Energy value must be positive.");
857  }
858  if (csv.real(i,1) > 0) {
859  log10value = std::log10(csv.real(i,1));
860  }
861  else {
863  csv.real(i,1), "Intensity value must be positive.");
864  }
865 
866  // Make sure that energies are increasing
867  if (csv.real(i,0) <= last_energy) {
869  csv.real(i,0), "Energy values must be monotonically increasing.");
870  }
871 
872  // Append log10 of node energy and value
873  m_lin_nodes.append(csv.real(i,0));
874  m_log_nodes.append(log10energy);
875  m_lin_values.push_back(csv.real(i,1));
876  m_log_values.push_back(log10value);
877 
878  // Store last energy for monotonically increasing check
879  last_energy = csv.real(i,0);
880 
881  } // endfor: looped over nodes
882 
883  // Set pre-computation cache
884  set_cache();
885 
886  } // endif: filename was not empty
887 
888  // Return
889  return;
890 }
891 
892 
893 /***********************************************************************//**
894  * @brief Set pre-computation cache
895  ***************************************************************************/
897 {
898  // Clear any existing values
899  m_prefactor.clear();
900  m_gamma.clear();
901  m_epivot.clear();
902  m_flux.clear();
903  m_eflux.clear();
904 
905  // Loop over all nodes-1
906  for (int i = 0; i < m_lin_nodes.size()-1; ++i) {
907 
908  // Get energies and function values
909  double emin = m_lin_nodes[i];
910  double emax = m_lin_nodes[i+1];
911  double fmin = m_lin_values[i];
912  double fmax = m_lin_values[i+1];
913 
914  // Compute pivot energy (MeV). We use here the geometric mean of
915  // the node boundaries.
916  double epivot = std::sqrt(emin*emax);
917 
918  // Compute spectral index
919  double gamma = std::log(fmin/fmax) / std::log(emin/emax);
920 
921  // Compute power law normalisation
922  double prefactor = fmin / std::pow(emin/epivot, gamma);
923 
924  // Compute photon flux between nodes
925  double flux = prefactor *
926  gammalib::plaw_photon_flux(emin, emax, epivot, gamma);
927 
928  // Compute energy flux between nodes
929  double eflux = prefactor *
930  gammalib::plaw_energy_flux(emin, emax, epivot, gamma);
931 
932  // Convert energy flux from MeV/cm2/s to erg/cm2/s
933  eflux *= gammalib::MeV2erg;
934 
935  // Push back values on pre-computation cache
936  m_prefactor.push_back(prefactor);
937  m_gamma.push_back(gamma);
938  m_epivot.push_back(epivot);
939  m_flux.push_back(flux);
940  m_eflux.push_back(eflux);
941 
942  } // endfor: looped over all nodes
943 
944  // Return
945  return;
946 }
947 
948 
949 /***********************************************************************//**
950  * @brief Set MC pre-computation cache
951  *
952  * @param[in] emin Minimum energy.
953  * @param[in] emax Maximum energy.
954  *
955  * This method sets up an array of indices and the cumulative distribution
956  * function needed for MC simulations.
957  ***************************************************************************/
959  const GEnergy& emax) const
960 {
961  // Check if we need to update the cache
962  if (emin != m_mc_emin || emax != m_mc_emax) {
963 
964  // Store new energy interval
965  m_mc_emin = emin;
966  m_mc_emax = emax;
967 
968  // Initialise cache
969  m_mc_cum.clear();
970  m_mc_min.clear();
971  m_mc_max.clear();
972  m_mc_exp.clear();
973 
974  // Get energy range in MeV
975  double e_min = emin.MeV();
976  double e_max = emax.MeV();
977 
978  // Continue only if e_max > e_min
979  if (e_max > e_min) {
980 
981  // Allocate flux
982  double flux;
983 
984  // Determine left node index for minimum energy
985  m_lin_nodes.set_value(e_min);
986  int inx_emin = m_lin_nodes.inx_left();
987 
988  // Determine left node index for maximum energy
989  m_lin_nodes.set_value(e_max);
990  int inx_emax = m_lin_nodes.inx_left();
991 
992  // If both energies are within the same node then just
993  // add this one node on the stack
994  if (inx_emin == inx_emax) {
995  flux = m_prefactor[inx_emin] *
997  e_max,
998  m_epivot[inx_emin],
999  m_gamma[inx_emin]);
1000  m_mc_cum.push_back(flux);
1001  m_mc_min.push_back(e_min);
1002  m_mc_max.push_back(e_max);
1003  m_mc_exp.push_back(m_gamma[inx_emin]);
1004  }
1005 
1006  // ... otherwise integrate over the nodes where emin and emax
1007  // resides and all the remaining nodes
1008  else {
1009 
1010  // If we are requested to extrapolate beyond first node,
1011  // the use the first nodes lower energy and upper energy
1012  // boundary for the initial integration.
1013  int i_start = (e_min < m_lin_nodes[0]) ? inx_emin : inx_emin+1;
1014 
1015  // Add emin to the node boundary
1016  flux = m_prefactor[inx_emin] *
1018  m_lin_nodes[i_start],
1019  m_epivot[inx_emin],
1020  m_gamma[inx_emin]);
1021  m_mc_cum.push_back(flux);
1022  m_mc_min.push_back(e_min);
1023  m_mc_max.push_back(m_lin_nodes[i_start]);
1024  m_mc_exp.push_back(m_gamma[inx_emin]);
1025 
1026  // Add all nodes between
1027  for (int i = i_start; i < inx_emax; ++i) {
1028  flux = m_flux[i];
1029  m_mc_cum.push_back(flux);
1030  m_mc_min.push_back(m_lin_nodes[i]);
1031  m_mc_max.push_back(m_lin_nodes[i+1]);
1032  m_mc_exp.push_back(m_gamma[i]);
1033  }
1034 
1035  // Add node boundary to emax
1036  flux = m_prefactor[inx_emax] *
1038  e_max,
1039  m_epivot[inx_emax],
1040  m_gamma[inx_emax]);
1041  m_mc_cum.push_back(flux);
1042  m_mc_min.push_back(m_lin_nodes[inx_emax]);
1043  m_mc_max.push_back(e_max);
1044  m_mc_exp.push_back(m_gamma[inx_emax]);
1045 
1046  } // endelse: emin and emax not between same nodes
1047 
1048  // Build cumulative distribution
1049  for (int i = 1; i < m_mc_cum.size(); ++i) {
1050  m_mc_cum[i] += m_mc_cum[i-1];
1051  }
1052  double norm = m_mc_cum[m_mc_cum.size()-1];
1053  if (norm > 0.0) {
1054  for (int i = 0; i < m_mc_cum.size(); ++i) {
1055  m_mc_cum[i] /= norm;
1056  }
1057  }
1058 
1059  // Set MC values
1060  for (int i = 0; i < m_mc_cum.size(); ++i) {
1061 
1062  // Compute exponent
1063  double exponent = m_mc_exp[i] + 1.0;
1064 
1065  // Exponent dependend computation
1066  if (std::abs(exponent) > 1.0e-11) {
1067 
1068  // If the exponent is too small then use lower energy
1069  // boundary
1070  if (exponent < -50.0) {
1071  m_mc_exp[i] = 0.0;
1072  m_mc_min[i] = std::log(m_mc_min[i]);
1073  m_mc_max[i] = m_mc_min[i];
1074  }
1075 
1076  // ... otherwise if exponent is too large then use
1077  // upper energy boundary
1078  else if (exponent > +50.0) {
1079  m_mc_exp[i] = 0.0;
1080  m_mc_min[i] = std::log(m_mc_max[i]);
1081  m_mc_max[i] = m_mc_min[i];
1082  }
1083 
1084  // ... otherwise use transformation formula
1085  else {
1086  m_mc_exp[i] = exponent;
1087  m_mc_min[i] = std::pow(m_mc_min[i], exponent);
1088  m_mc_max[i] = std::pow(m_mc_max[i], exponent);
1089  }
1090  }
1091  else {
1092  m_mc_exp[i] = 0.0;
1093  m_mc_min[i] = std::log(m_mc_min[i]);
1094  m_mc_max[i] = std::log(m_mc_max[i]);
1095  }
1096 
1097  } // endfor: set MC values
1098 
1099  } // endif: e_max > e_min
1100 
1101  } // endif: Update was required
1102 
1103 
1104  // Return
1105  return;
1106 }
std::vector< double > m_gamma
Power-law indices.
void mc_update(const GEnergy &emin, const GEnergy &emax) const
Set MC pre-computation cache.
GEnergy m_mc_emax
Maximum energy.
int size(void) const
Return number of nodes in node array.
Definition: GNodeArray.hpp:188
virtual double flux(const GEnergy &emin, const GEnergy &emax) const
Returns model photon flux between emin, emax
std::vector< double > m_lin_values
Function values at nodes.
const double & factor_gradient(void) const
Return parameter gradient factor.
double norm(const GVector &vector)
Computes vector norm.
Definition: GVector.cpp:821
const std::string & name(void) const
Return parameter name.
Random number generator class definition.
Abstract spectral model base class.
double gradient(void) const
Return parameter gradient.
#define G_WRITE
std::vector< double > m_mc_max
Upper boundary for MC.
int size(void) const
Return number of parameters.
GVector abs(const GVector &vector)
Computes absolute of vector elements.
Definition: GVector.cpp:1163
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.
bool is_empty(void) const
Signal if filename is empty.
Definition: GFilename.hpp:160
std::vector< double > m_mc_exp
Exponent for MC.
XML element node class.
Definition: GXmlElement.hpp:47
virtual void read(const GXmlElement &xml)
Read model from XML element.
Spectral model registry class definition.
void clear(void)
Clear node array.
Definition: GNodeArray.cpp:235
Random number generator class.
Definition: GRan.hpp:44
Comma-separated values table class.
Definition: GCsv.hpp:57
#define G_MC
double MeV(void) const
Return energy in MeV.
Definition: GEnergy.cpp:321
virtual int elements(void) const
Return number of GXMLElement children of node.
Definition: GXmlNode.cpp:580
Time class.
Definition: GTime.hpp:54
void load_nodes(const GFilename &filename)
Load nodes from file.
Gammalib tools definition.
void free_members(void)
Delete class members.
void set_value(const double &value) const
Set indices and weighting factors for interpolation.
Definition: GNodeArray.cpp:580
GFilename xml_file_reduce(const GXmlElement &xml, const std::string &filename)
Reduce file name provided for writing as XML attribute.
Definition: GTools.cpp:1713
const GFilename & filename(void) const
Return file name.
double plaw_energy_flux(const double &emin, const double &emax, const double &epivot, const double &gamma)
Compute energy flux between two energies for a power law.
Definition: GTools.cpp:1127
bool is_free(void) const
Signal if parameter is free.
std::vector< double > m_log_values
log10(Function) values at nodes
virtual void write(GXmlElement &xml) const
Write model into XML element.
double log10MeV(void) const
Return log10 of energy in MeV.
Definition: GEnergy.cpp:423
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.
GFilename m_filename
Name of file function.
bool is_infinite(const double &x)
Signal if argument is infinite.
Definition: GTools.hpp:167
virtual ~GModelSpectralFunc(void)
Destructor.
const int & ncols(void) const
Return number of columns.
Definition: GCsv.hpp:137
std::vector< double > m_eflux
Energy fluxes.
GVector sqrt(const GVector &vector)
Computes square root of vector elements.
Definition: GVector.cpp:1268
bool has_grad(void) const
Signal if parameter gradient is computed analytically.
const GXmlAttribute * attribute(const int &index) const
Return attribute.
std::vector< double > m_mc_min
Lower boundary for MC.
GEnergy m_mc_emin
Minimum energy.
void free(void)
Free a parameter.
virtual void clear(void)
Clear file function.
Filename class.
Definition: GFilename.hpp:62
void set_cache(void) const
Set pre-computation cache.
GNodeArray m_log_nodes
lof10(Energy) nodes of function
Spectral function model class definition.
GModelPar m_norm
Normalization factor.
GVector log(const GVector &vector)
Computes natural logarithm of vector elements.
Definition: GVector.cpp:1184
double uniform(void)
Returns random double precision floating value in range 0 to 1.
Definition: GRan.cpp:242
void clear(void)
Clear parameter.
virtual double eflux(const GEnergy &emin, const GEnergy &emax) const
Returns model energy flux between emin, emax
double interpolate(const double &value, const std::vector< double > &vector) const
Interpolate value.
Definition: GNodeArray.cpp:536
#define G_LOAD_NODES
GChatter
Definition: GTypemaps.hpp:33
const int & inx_left(void) const
Returns left node index.
Definition: GNodeArray.hpp:231
const int & nrows(void) const
Return number of rows.
Definition: GCsv.hpp:149
void copy_members(const GModelSpectralFunc &model)
Copy class members.
Interface definition for the spectral model registry class.
std::vector< double > m_epivot
Power-law pivot energies.
double real(const int &row, const int &col) const
Get double precision value.
Definition: GCsv.cpp:318
Spectral function model class.
std::string url(void) const
Return Uniform Resource Locator (URL)
Definition: GFilename.hpp:189
virtual GXmlElement * element(const int &index)
Return pointer to GXMLElement child.
Definition: GXmlNode.cpp:634
void clear(void)
Clear file name.
Definition: GFilename.cpp:188
void init_members(void)
Initialise class members.
virtual GEnergy mc(const GEnergy &emin, const GEnergy &emax, const GTime &time, GRan &ran) const
Returns MC energy between [emin, emax].
const double MeV2erg
Definition: GTools.hpp:45
std::vector< double > m_mc_cum
Cumulative distribution.
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
virtual GModelSpectralFunc * clone(void) const
Clone file function.
double value(void) const
Return parameter value.
GVector pow(const GVector &vector, const double &power)
Computes tanh of vector elements.
Definition: GVector.cpp:1332
std::vector< double > m_prefactor
Power-law normalisations.
Exception handler interface definition.
virtual std::string print(const GChatter &chatter=NORMAL) const
Print file function information.
void init_members(void)
Initialise class members.
GVector exp(const GVector &vector)
Computes exponential of vector elements.
Definition: GVector.cpp:1142
#define G_READ
GModelSpectralFunc(void)
Void constructor.
GNodeArray m_lin_nodes
Energy nodes of function.
virtual GXmlNode * append(const GXmlNode &node)
Append XML child node.
Definition: GXmlNode.cpp:284
std::string parformat(const std::string &s, const int &indent=0)
Convert string in parameter format.
Definition: GTools.cpp:1022
void append(const double &node)
Append one node to array.
Definition: GNodeArray.cpp:351
const GModelSpectralFunc g_spectral_func_seed
virtual GModelSpectralFunc & operator=(const GModelSpectralFunc &model)
Assignment operator.
double norm(void) const
Return normalization factor.
void free_members(void)
Delete class members.
void clear(void)
Clear instance.
Definition: GEnergy.cpp:261
virtual double eval(const GEnergy &srcEng, const GTime &srcTime=GTime(), const bool &gradients=false) const
Evaluate function.
Comma-separated values table class definition.
GVector log10(const GVector &vector)
Computes base10 logarithm of vector elements.
Definition: GVector.cpp:1205
Class that handles energies in a unit independent way.
Definition: GEnergy.hpp:48
double plaw_photon_flux(const double &emin, const double &emax, const double &epivot, const double &gamma)
Compute photon flux between two energies for a power law.
Definition: GTools.cpp:1079
GFilename xml_file_expand(const GXmlElement &xml, const std::string &filename)
Expand file name provided as XML attribute for loading.
Definition: GTools.cpp:1683
std::vector< double > m_flux
Photon fluxes.
std::string str(const unsigned short int &value)
Convert unsigned short integer value into string.
Definition: GTools.cpp:413