GammaLib  2.0.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} F(E)
244  * \f]
245  *
246  * where
247  * - \f${\tt m\_norm}\f$ is the normalization factor and
248  * - \f${F(E)}\f$ is the spectral function.
249  *
250  * If the @p gradients flag is true the method will also compute the
251  * partial derivatives of the model with respect to the parameters using
252  *
253  * \f[
254  * \frac{\delta S_{\rm E}(E | t)}{\delta {\tt m\_norm}} =
255  * \frac{S_{\rm E}(E | t)}{{\tt m\_norm}}
256  * \f]
257  ***************************************************************************/
258 double GModelSpectralFunc::eval(const GEnergy& srcEng,
259  const GTime& srcTime,
260  const bool& gradients) const
261 {
262  // Interpolate function. This is done in log10-log10 space, but the
263  // linear value is returned.
264  double arg = m_log_nodes.interpolate(srcEng.log10MeV(), m_log_values);
265  double func = std::pow(10.0, arg);
266 
267  // Compute function value
268  double value = m_norm.value() * func;
269 
270  // Optionally compute gradients
271  if (gradients) {
272 
273  // Compute partial derivatives of the parameter values
274  double g_norm = (m_norm.is_free()) ? m_norm.scale() * func : 0.0;
275 
276  // Set gradients
277  m_norm.factor_gradient(g_norm);
278 
279  } // endif: gradient computation was requested
280 
281  // Compile option: Check for NaN/Inf
282  #if defined(G_NAN_CHECK)
283  if (gammalib::is_notanumber(value) || gammalib::is_infinite(value)) {
284  std::cout << "*** ERROR: GModelSpectralFunc::eval";
285  std::cout << "(srcEng=" << srcEng;
286  std::cout << ", srcTime=" << srcTime << "):";
287  std::cout << " NaN/Inf encountered";
288  std::cout << " (value=" << value;
289  std::cout << ", norm=" << norm();
290  std::cout << ", func=" << func;
291  std::cout << ")" << std::endl;
292  }
293  #endif
294 
295  // Return
296  return value;
297 }
298 
299 
300 /***********************************************************************//**
301  * @brief Returns model photon flux between [emin, emax] (units: 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  ***************************************************************************/
317 double GModelSpectralFunc::flux(const GEnergy& emin, const GEnergy& emax) const
318 {
319  // Initialise flux
320  double flux = 0.0;
321 
322  // Compute only if integration range is valid
323  if (emin < emax) {
324 
325  // Get energy range in MeV
326  double e_min = emin.MeV();
327  double e_max = emax.MeV();
328 
329  // Determine left node index for minimum energy
330  m_lin_nodes.set_value(e_min);
331  int inx_emin = m_lin_nodes.inx_left();
332 
333  // Determine left node index for maximum energy
334  m_lin_nodes.set_value(e_max);
335  int inx_emax = m_lin_nodes.inx_left();
336 
337  // If both energies are within the same nodes then simply
338  // integrate over the energy interval using the appropriate power
339  // law parameters
340  if (inx_emin == inx_emax) {
341  flux = m_prefactor[inx_emin] *
343  e_max,
344  m_epivot[inx_emin],
345  m_gamma[inx_emin]);
346  }
347 
348  // ... otherwise integrate over the nodes where emin and emax
349  // resides and all the remaining nodes
350  else {
351 
352  // If we are requested to extrapolate beyond first node,
353  // the use the first nodes lower energy and upper energy
354  // boundary for the initial integration.
355  int i_start = (e_min < m_lin_nodes[0]) ? inx_emin : inx_emin+1;
356 
357  // Integrate from emin to the node boundary
358  flux = m_prefactor[inx_emin] *
360  m_lin_nodes[i_start],
361  m_epivot[inx_emin],
362  m_gamma[inx_emin]);
363 
364  // Integrate over all nodes between
365  for (int i = i_start; i < inx_emax; ++i) {
366  flux += m_flux[i];
367  }
368 
369  // Integrate from node boundary to emax
370  flux += m_prefactor[inx_emax] *
372  e_max,
373  m_epivot[inx_emax],
374  m_gamma[inx_emax]);
375 
376  } // endelse: emin and emax not between same nodes
377 
378  // Multiply flux by normalisation factor
379  flux *= norm();
380 
381  } // endif: integration range was valid
382 
383  // Return
384  return flux;
385 }
386 
387 
388 /***********************************************************************//**
389  * @brief Returns model energy flux between [emin, emax] (units: erg/cm2/s)
390  *
391  * @param[in] emin Minimum photon energy.
392  * @param[in] emax Maximum photon energy.
393  * @return Energy flux (erg/cm2/s).
394  *
395  * Computes
396  *
397  * \f[
398  * \int_{\tt emin}^{\tt emax} S_{\rm E}(E | t) E \, dE
399  * \f]
400  *
401  * where
402  * - [@p emin, @p emax] is an energy interval, and
403  * - \f$S_{\rm E}(E | t)\f$ is the spectral model (ph/cm2/s/MeV).
404  ***************************************************************************/
405 double GModelSpectralFunc::eflux(const GEnergy& emin, const GEnergy& emax) const
406 {
407  // Initialise flux
408  double eflux = 0.0;
409 
410  // Compute only if integration range is valid
411  if (emin < emax) {
412 
413  // Get energy range in MeV
414  double e_min = emin.MeV();
415  double e_max = emax.MeV();
416 
417  // Determine left node index for minimum energy
418  m_lin_nodes.set_value(e_min);
419  int inx_emin = m_lin_nodes.inx_left();
420 
421  // Determine left node index for maximum energy
422  m_lin_nodes.set_value(e_max);
423  int inx_emax = m_lin_nodes.inx_left();
424 
425  // If both energies are within the same nodes then simply
426  // integrate over the energy interval using the appropriate power
427  // law parameters
428  if (inx_emin == inx_emax) {
429  eflux = m_prefactor[inx_emin] *
431  e_max,
432  m_epivot[inx_emin],
433  m_gamma[inx_emin]) *
435  }
436 
437  // ... otherwise integrate over the nodes where emin and emax
438  // resides and all the remaining nodes
439  else {
440 
441  // If we are requested to extrapolate beyond first node,
442  // the use the first nodes lower energy and upper energy
443  // boundary for the initial integration.
444  int i_start = (e_min < m_lin_nodes[0]) ? inx_emin : inx_emin+1;
445 
446  // Integrate from emin to the node boundary
447  eflux = m_prefactor[inx_emin] *
449  m_lin_nodes[i_start],
450  m_epivot[inx_emin],
451  m_gamma[inx_emin]) *
453 
454  // Integrate over all nodes between
455  for (int i = i_start; i < inx_emax; ++i) {
456  eflux += m_eflux[i];
457  }
458 
459  // Integrate from node boundary to emax
460  eflux += m_prefactor[inx_emax] *
462  e_max,
463  m_epivot[inx_emax],
464  m_gamma[inx_emax]) *
466 
467  } // endelse: emin and emax not between same nodes
468 
469  // Multiply flux by normalisation factor
470  eflux *= norm();
471 
472  } // endif: integration range was valid
473 
474  // Return flux
475  return eflux;
476 }
477 
478 
479 /***********************************************************************//**
480  * @brief Returns MC energy between [emin, emax]
481  *
482  * @param[in] emin Minimum photon energy.
483  * @param[in] emax Maximum photon energy.
484  * @param[in] time True photon arrival time.
485  * @param[in,out] ran Random number generator.
486  * @return Energy.
487  *
488  * @exception GException::erange_invalid
489  * Energy range is invalid (emin < emax required).
490  *
491  * Returns Monte Carlo energy by randomly drawing from a broken power law
492  * defined by the file function.
493  ***************************************************************************/
495  const GEnergy& emax,
496  const GTime& time,
497  GRan& ran) const
498 {
499  // Throw an exception if energy range is invalid
500  if (emin >= emax) {
501  throw GException::erange_invalid(G_MC, emin.MeV(), emax.MeV(),
502  "Minimum energy < maximum energy required.");
503  }
504 
505  // Allocate energy
506  GEnergy energy;
507 
508  // Update cache
509  mc_update(emin, emax);
510 
511  // Determine in which bin we reside
512  int inx = 0;
513  if (m_mc_cum.size() > 1) {
514  double u = ran.uniform();
515  for (inx = m_mc_cum.size()-1; inx > 0; --inx) {
516  if (m_mc_cum[inx-1] <= u) {
517  break;
518  }
519  }
520  }
521 
522  // Get random energy for specific bin
523  if (m_mc_exp[inx] != 0.0) {
524  double e_min = m_mc_min[inx];
525  double e_max = m_mc_max[inx];
526  double u = ran.uniform();
527  double eng = (u > 0.0)
528  ? std::exp(std::log(u * (e_max - e_min) + e_min) / m_mc_exp[inx])
529  : 0.0;
530  energy.MeV(eng);
531  }
532  else {
533  double e_min = m_mc_min[inx];
534  double e_max = m_mc_max[inx];
535  double u = ran.uniform();
536  double eng = std::exp(u * (e_max - e_min) + e_min);
537  energy.MeV(eng);
538  }
539 
540  // Return energy
541  return energy;
542 }
543 
544 
545 /***********************************************************************//**
546  * @brief Read model from XML element
547  *
548  * @param[in] xml XML element containing power law model information.
549  *
550  * @exception GException::model_invalid_parnum
551  * Invalid number of model parameters found in XML element.
552  * @exception GException::model_invalid_parnames
553  * Invalid model parameter name found in XML element.
554  *
555  * Reads the spectral information from an XML element. The format of the XML
556  * elements is
557  *
558  * <spectrum type="FileFunction" file="..">
559  * <parameter name="Normalization" scale=".." value=".." min=".." max=".." free=".."/>
560  * </spectrum>
561  ***************************************************************************/
563 {
564  // Verify that XML element has exactly 1 parameter
565  if (xml.elements() != 1 || xml.elements("parameter") != 1) {
567  "Spectral function requires exactly 1 parameter.");
568  }
569 
570  // Get parameter element
571  const GXmlElement* par = xml.element("parameter", 0);
572 
573  // Get value
574  if (par->attribute("name") == "Normalization") {
575  m_norm.read(*par);
576  }
577  else {
579  "Require \"Normalization\" parameter.");
580  }
581 
582  // Load nodes from file
584 
585  // Return
586  return;
587 }
588 
589 
590 /***********************************************************************//**
591  * @brief Write model into XML element
592  *
593  * @param[in] xml XML element into which model information is written.
594  *
595  * @exception GException::model_invalid_spectral
596  * Existing XML element is not of type "FileFunction"
597  * @exception GException::model_invalid_parnum
598  * Invalid number of model parameters or nodes found in XML element.
599  * @exception GException::model_invalid_parnames
600  * Invalid model parameter names found in XML element.
601  *
602  * Writes the spectral information into an XML element. The format of the XML
603  * element is
604  *
605  * <spectrum type="FileFunction" file="..">
606  * <parameter name="Normalization" scale=".." value=".." min=".." max=".." free=".."/>
607  * </spectrum>
608  *
609  * Note that the function nodes will not be written since they will not be
610  * altered by any method.
611  ***************************************************************************/
613 {
614  // Set model type
615  if (xml.attribute("type") == "") {
616  xml.attribute("type", "FileFunction");
617  }
618 
619  // Verify model type
620  if (xml.attribute("type") != "FileFunction") {
622  "Spectral model is not of type \"FileFunction\".");
623  }
624 
625  // If XML element has 0 nodes then append 1 parameter node
626  if (xml.elements() == 0) {
627  xml.append(GXmlElement("parameter name=\"Normalization\""));
628  }
629 
630  // Verify that XML element has exactly 1 parameter
631  if (xml.elements() != 1 || xml.elements("parameter") != 1) {
633  "Spectral function requires exactly 1 parameter.");
634  }
635 
636  // Get parameter element
637  GXmlElement* par = xml.element("parameter", 0);
638 
639  // Set parameter
640  if (par->attribute("name") == "Normalization") {
641  m_norm.write(*par);
642  }
643  else {
645  "Require \"Normalization\" parameter.");
646  }
647 
648  // Set file attribute
650 
651  // Return
652  return;
653 }
654 
655 
656 /***********************************************************************//**
657  * @brief Print file function information
658  *
659  * @param[in] chatter Chattiness (defaults to NORMAL).
660  * @return String containing file function information.
661  ***************************************************************************/
662 std::string GModelSpectralFunc::print(const GChatter& chatter) const
663 {
664  // Initialise result string
665  std::string result;
666 
667  // Continue only if chatter is not silent
668  if (chatter != SILENT) {
669 
670  // Append header
671  result.append("=== GModelSpectralFunc ===");
672 
673  // Append information
674  result.append("\n"+gammalib::parformat("Function file"));
675  result.append(m_filename.url());
676  result.append("\n"+gammalib::parformat("Number of nodes"));
677  result.append(gammalib::str(m_lin_nodes.size()));
678  result.append("\n"+gammalib::parformat("Number of parameters"));
679  result.append(gammalib::str(size()));
680  for (int i = 0; i < size(); ++i) {
681  result.append("\n"+m_pars[i]->print(chatter));
682  }
683 
684  // Append node information
685  for (int i = 0; i < m_prefactor.size(); ++i) {
686  result.append("\n"+gammalib::parformat("Node "+gammalib::str(i+1)));
687  result.append("Epivot="+gammalib::str(m_epivot[i]));
688  result.append(" Prefactor="+gammalib::str(m_prefactor[i]));
689  result.append(" Gamma="+gammalib::str(m_gamma[i]));
690  result.append(" Flux="+gammalib::str(m_flux[i]));
691  result.append(" EFlux="+gammalib::str(m_eflux[i]));
692  }
693 
694  } // endif: chatter was not silent
695 
696  // Return result
697  return result;
698 }
699 
700 
701 /*==========================================================================
702  = =
703  = Private methods =
704  = =
705  ==========================================================================*/
706 
707 /***********************************************************************//**
708  * @brief Initialise class members
709  ***************************************************************************/
711 {
712  // Initialise powerlaw normalisation
713  m_norm.clear();
714  m_norm.name("Normalization");
715  m_norm.scale(1.0);
716  m_norm.value(1.0);
717  m_norm.range(0.0,1000.0);
718  m_norm.free();
719  m_norm.gradient(0.0);
720  m_norm.has_grad(true);
721 
722  // Set parameter pointer(s)
723  m_pars.clear();
724  m_pars.push_back(&m_norm);
725 
726  // Initialise members
727  m_lin_nodes.clear();
728  m_log_nodes.clear();
729  m_lin_values.clear();
730  m_log_values.clear();
731  m_filename.clear();
732 
733  // Initialise flux cache
734  m_prefactor.clear();
735  m_gamma.clear();
736  m_epivot.clear();
737  m_flux.clear();
738  m_eflux.clear();
739 
740  // Initialise MC cache
741  m_mc_emin.clear();
742  m_mc_emax.clear();
743  m_mc_cum.clear();
744  m_mc_min.clear();
745  m_mc_max.clear();
746  m_mc_exp.clear();
747 
748  // Return
749  return;
750 }
751 
752 
753 /***********************************************************************//**
754  * @brief Copy class members
755  *
756  * @param[in] model Spectral function model.
757  ***************************************************************************/
759 {
760  // Copy members
761  m_norm = model.m_norm;
762  m_lin_nodes = model.m_lin_nodes;
763  m_log_nodes = model.m_log_nodes;
764  m_lin_values = model.m_lin_values;
765  m_log_values = model.m_log_values;
766  m_filename = model.m_filename;
767 
768  // Copy flux cache
769  m_prefactor = model.m_prefactor;
770  m_gamma = model.m_gamma;
771  m_epivot = model.m_epivot;
772  m_flux = model.m_flux;
773  m_eflux = model.m_eflux;
774 
775  // Copy MC cache
776  m_mc_emin = model.m_mc_emin;
777  m_mc_emax = model.m_mc_emax;
778  m_mc_cum = model.m_mc_cum;
779  m_mc_min = model.m_mc_min;
780  m_mc_max = model.m_mc_max;
781  m_mc_exp = model.m_mc_exp;
782 
783  // Set parameter pointer(s)
784  m_pars.clear();
785  m_pars.push_back(&m_norm);
786 
787  // Return
788  return;
789 }
790 
791 
792 /***********************************************************************//**
793  * @brief Delete class members
794  ***************************************************************************/
796 {
797  // Return
798  return;
799 }
800 
801 
802 /***********************************************************************//**
803  * @brief Load nodes from file
804  *
805  * @param[in] filename File name.
806  *
807  * @exception GException::file_function_data
808  * File contains less than 2 nodes
809  * @exception GException::file_function_columns
810  * File contains less than 2 columns
811  * @exception GException::file_function_value
812  * File contains invalid value
813  *
814  * The file function is stored as a column separated value table (CSV) in an
815  * ASCII file with (at least) 2 columns. The first column specifies the
816  * energy in MeV while the second column specifies the intensity at this
817  * energy in units of ph/cm2/s/MeV.
818  * The node energies and values will be stored both linearly and as log10.
819  * The log10 storing requires that node energies and node values are
820  * positive. Also, at least 2 nodes and 2 columns are required in the file
821  * function.
822  ***************************************************************************/
824 {
825  // Clear nodes and values
826  m_lin_nodes.clear();
827  m_log_nodes.clear();
828  m_lin_values.clear();
829  m_log_values.clear();
830 
831  // Set filename
833 
834  // Continue only if filename is not empty
835  if (!filename.is_empty()) {
836 
837  // Load file
838  GCsv csv = GCsv(filename.url());
839 
840  // Check if there are at least 2 nodes
841  if (csv.nrows() < 2) {
843  csv.nrows());
844  }
845 
846  // Check if there are at least 2 columns
847  if (csv.ncols() < 2) {
849  csv.ncols());
850  }
851 
852  // Setup nodes
853  double last_energy = 0.0;
854  for (int i = 0; i < csv.nrows(); ++i) {
855 
856  // Get log10 of node energy and value. Make sure they are valid.
857  double log10energy;
858  double log10value;
859  if (csv.real(i,0) > 0) {
860  log10energy = std::log10(csv.real(i,0));
861  }
862  else {
864  csv.real(i,0), "Energy value must be positive.");
865  }
866  if (csv.real(i,1) > 0) {
867  log10value = std::log10(csv.real(i,1));
868  }
869  else {
871  csv.real(i,1), "Intensity value must be positive.");
872  }
873 
874  // Make sure that energies are increasing
875  if (csv.real(i,0) <= last_energy) {
877  csv.real(i,0), "Energy values must be monotonically increasing.");
878  }
879 
880  // Append log10 of node energy and value
881  m_lin_nodes.append(csv.real(i,0));
882  m_log_nodes.append(log10energy);
883  m_lin_values.push_back(csv.real(i,1));
884  m_log_values.push_back(log10value);
885 
886  // Store last energy for monotonically increasing check
887  last_energy = csv.real(i,0);
888 
889  } // endfor: looped over nodes
890 
891  // Set pre-computation cache
892  set_cache();
893 
894  } // endif: filename was not empty
895 
896  // Return
897  return;
898 }
899 
900 
901 /***********************************************************************//**
902  * @brief Set pre-computation cache
903  ***************************************************************************/
905 {
906  // Clear any existing values
907  m_prefactor.clear();
908  m_gamma.clear();
909  m_epivot.clear();
910  m_flux.clear();
911  m_eflux.clear();
912 
913  // Loop over all nodes-1
914  for (int i = 0; i < m_lin_nodes.size()-1; ++i) {
915 
916  // Get energies and function values
917  double emin = m_lin_nodes[i];
918  double emax = m_lin_nodes[i+1];
919  double fmin = m_lin_values[i];
920  double fmax = m_lin_values[i+1];
921 
922  // Compute pivot energy (MeV). We use here the geometric mean of
923  // the node boundaries.
924  double epivot = std::sqrt(emin*emax);
925 
926  // Compute spectral index
927  double gamma = std::log(fmin/fmax) / std::log(emin/emax);
928 
929  // Compute power law normalisation
930  double prefactor = fmin / std::pow(emin/epivot, gamma);
931 
932  // Compute photon flux between nodes
933  double flux = prefactor *
934  gammalib::plaw_photon_flux(emin, emax, epivot, gamma);
935 
936  // Compute energy flux between nodes
937  double eflux = prefactor *
938  gammalib::plaw_energy_flux(emin, emax, epivot, gamma);
939 
940  // Convert energy flux from MeV/cm2/s to erg/cm2/s
941  eflux *= gammalib::MeV2erg;
942 
943  // Push back values on pre-computation cache
944  m_prefactor.push_back(prefactor);
945  m_gamma.push_back(gamma);
946  m_epivot.push_back(epivot);
947  m_flux.push_back(flux);
948  m_eflux.push_back(eflux);
949 
950  } // endfor: looped over all nodes
951 
952  // Return
953  return;
954 }
955 
956 
957 /***********************************************************************//**
958  * @brief Set MC pre-computation cache
959  *
960  * @param[in] emin Minimum energy.
961  * @param[in] emax Maximum energy.
962  *
963  * This method sets up an array of indices and the cumulative distribution
964  * function needed for MC simulations.
965  ***************************************************************************/
967  const GEnergy& emax) const
968 {
969  // Check if we need to update the cache
970  if (emin != m_mc_emin || emax != m_mc_emax) {
971 
972  // Store new energy interval
973  m_mc_emin = emin;
974  m_mc_emax = emax;
975 
976  // Initialise cache
977  m_mc_cum.clear();
978  m_mc_min.clear();
979  m_mc_max.clear();
980  m_mc_exp.clear();
981 
982  // Get energy range in MeV
983  double e_min = emin.MeV();
984  double e_max = emax.MeV();
985 
986  // Continue only if e_max > e_min
987  if (e_max > e_min) {
988 
989  // Allocate flux
990  double flux;
991 
992  // Determine left node index for minimum energy
993  m_lin_nodes.set_value(e_min);
994  int inx_emin = m_lin_nodes.inx_left();
995 
996  // Determine left node index for maximum energy
997  m_lin_nodes.set_value(e_max);
998  int inx_emax = m_lin_nodes.inx_left();
999 
1000  // If both energies are within the same node then just
1001  // add this one node on the stack
1002  if (inx_emin == inx_emax) {
1003  flux = m_prefactor[inx_emin] *
1005  e_max,
1006  m_epivot[inx_emin],
1007  m_gamma[inx_emin]);
1008  m_mc_cum.push_back(flux);
1009  m_mc_min.push_back(e_min);
1010  m_mc_max.push_back(e_max);
1011  m_mc_exp.push_back(m_gamma[inx_emin]);
1012  }
1013 
1014  // ... otherwise integrate over the nodes where emin and emax
1015  // resides and all the remaining nodes
1016  else {
1017 
1018  // If we are requested to extrapolate beyond first node,
1019  // the use the first nodes lower energy and upper energy
1020  // boundary for the initial integration.
1021  int i_start = (e_min < m_lin_nodes[0]) ? inx_emin : inx_emin+1;
1022 
1023  // Add emin to the node boundary
1024  flux = m_prefactor[inx_emin] *
1026  m_lin_nodes[i_start],
1027  m_epivot[inx_emin],
1028  m_gamma[inx_emin]);
1029  m_mc_cum.push_back(flux);
1030  m_mc_min.push_back(e_min);
1031  m_mc_max.push_back(m_lin_nodes[i_start]);
1032  m_mc_exp.push_back(m_gamma[inx_emin]);
1033 
1034  // Add all nodes between
1035  for (int i = i_start; i < inx_emax; ++i) {
1036  flux = m_flux[i];
1037  m_mc_cum.push_back(flux);
1038  m_mc_min.push_back(m_lin_nodes[i]);
1039  m_mc_max.push_back(m_lin_nodes[i+1]);
1040  m_mc_exp.push_back(m_gamma[i]);
1041  }
1042 
1043  // Add node boundary to emax
1044  flux = m_prefactor[inx_emax] *
1046  e_max,
1047  m_epivot[inx_emax],
1048  m_gamma[inx_emax]);
1049  m_mc_cum.push_back(flux);
1050  m_mc_min.push_back(m_lin_nodes[inx_emax]);
1051  m_mc_max.push_back(e_max);
1052  m_mc_exp.push_back(m_gamma[inx_emax]);
1053 
1054  } // endelse: emin and emax not between same nodes
1055 
1056  // Build cumulative distribution
1057  for (int i = 1; i < m_mc_cum.size(); ++i) {
1058  m_mc_cum[i] += m_mc_cum[i-1];
1059  }
1060  double norm = m_mc_cum[m_mc_cum.size()-1];
1061  if (norm > 0.0) {
1062  for (int i = 0; i < m_mc_cum.size(); ++i) {
1063  m_mc_cum[i] /= norm;
1064  }
1065  }
1066 
1067  // Set MC values
1068  for (int i = 0; i < m_mc_cum.size(); ++i) {
1069 
1070  // Compute exponent
1071  double exponent = m_mc_exp[i] + 1.0;
1072 
1073  // Exponent dependend computation
1074  if (std::abs(exponent) > 1.0e-11) {
1075 
1076  // If the exponent is too small then use lower energy
1077  // boundary
1078  if (exponent < -50.0) {
1079  m_mc_exp[i] = 0.0;
1080  m_mc_min[i] = std::log(m_mc_min[i]);
1081  m_mc_max[i] = m_mc_min[i];
1082  }
1083 
1084  // ... otherwise if exponent is too large then use
1085  // upper energy boundary
1086  else if (exponent > +50.0) {
1087  m_mc_exp[i] = 0.0;
1088  m_mc_min[i] = std::log(m_mc_max[i]);
1089  m_mc_max[i] = m_mc_min[i];
1090  }
1091 
1092  // ... otherwise use transformation formula
1093  else {
1094  m_mc_exp[i] = exponent;
1095  m_mc_min[i] = std::pow(m_mc_min[i], exponent);
1096  m_mc_max[i] = std::pow(m_mc_max[i], exponent);
1097  }
1098  }
1099  else {
1100  m_mc_exp[i] = 0.0;
1101  m_mc_min[i] = std::log(m_mc_min[i]);
1102  m_mc_max[i] = std::log(m_mc_max[i]);
1103  }
1104 
1105  } // endfor: set MC values
1106 
1107  } // endif: e_max > e_min
1108 
1109  } // endif: Update was required
1110 
1111  // Return
1112  return;
1113 }
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:192
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 factor gradient.
double norm(const GVector &vector)
Computes vector norm.
Definition: GVector.cpp:828
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:1170
virtual GModelSpectral & operator=(const GModelSpectral &model)
Assignment operator.
void write(GXmlElement &xml) const
Set or update parameter attributes in XML element.
Definition: GModelPar.cpp:351
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:1733
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:1138
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:1275
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
log10(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:1191
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:235
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:1339
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:1149
#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:1033
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:1212
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:1090
GFilename xml_file_expand(const GXmlElement &xml, const std::string &filename)
Expand file name provided as XML attribute for loading.
Definition: GTools.cpp:1703
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:415