GammaLib  2.0.0
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
GModelSky.cpp
Go to the documentation of this file.
1 /***************************************************************************
2  * GModelSky.cpp - Sky model class *
3  * ----------------------------------------------------------------------- *
4  * copyright (C) 2011-2021 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 GModelSky.cpp
23  * @brief Sky model class implementation
24  * @author Juergen Knoedlseder
25  */
26 
27 /* __ Includes ___________________________________________________________ */
28 #ifdef HAVE_CONFIG_H
29 #include <config.h>
30 #endif
31 #include "GTools.hpp"
32 #include "GException.hpp"
33 #include "GIntegral.hpp"
34 #include "GEnergies.hpp"
35 #include "GModelRegistry.hpp"
36 #include "GModelSky.hpp"
39 #include "GModelSpatialRadial.hpp"
45 #include "GModelTemporalConst.hpp"
46 #include "GSource.hpp"
47 #include "GResponse.hpp"
48 
49 /* __ Globals ____________________________________________________________ */
50 const GModelSky g_pointsource_seed("PointSource");
51 const GModelSky g_extendedsource_seed("ExtendedSource");
52 const GModelSky g_diffusesource_seed("DiffuseSource");
53 const GModelSky g_compositesource_seed("CompositeSource");
54 const GModelRegistry g_pointsource_registry(&g_pointsource_seed);
55 const GModelRegistry g_extendedsource_registry(&g_extendedsource_seed);
56 const GModelRegistry g_diffusesource_registry(&g_diffusesource_seed);
57 const GModelRegistry g_compositesource_registry(&g_compositesource_seed);
58 
59 /* __ Method name definitions ____________________________________________ */
60 #define G_NPRED "GModelSky::npred(GEnergy&, GTime&, GObservation&)"
61 #define G_XML_SPATIAL "GModelSky::xml_spatial(GXmlElement&)"
62 #define G_XML_SPECTRAL "GModelSky::xml_spectral(GXmlElement&)"
63 #define G_XML_TEMPORAL "GModelSky::xml_temporal(GXmlElement&)"
64 #define G_INTEGRATE_TIME "GModelSky::integrate_time(GEvent&, GObservation&,"\
65  " bool)"
66 #define G_INTEGRATE_ENERGY "GModelSky::integrate_energy(GEvent&, GTime&," \
67  " GObservation&, bool)"
68 #define G_INTEGRATE_DIR "GModelSky::integrate_dir(GEvent&, GEnergy&, GTime&,"\
69  " GObservation, bool)"
70 
71 /* __ Macros _____________________________________________________________ */
72 
73 /* __ Coding definitions _________________________________________________ */
74 
75 /* __ Debug definitions __________________________________________________ */
76 //#define G_DUMP_MC //!< Dump MC information
77 //#define G_DUMP_MC_DETAIL //!< Dump detailed MC information
78 
79 
80 /*==========================================================================
81  = =
82  = Constructors/destructors =
83  = =
84  ==========================================================================*/
85 
86 /***********************************************************************//**
87  * @brief Void constructor
88  *
89  * Construct an empty sky model. An empty sky model has no type.
90  ***************************************************************************/
92 {
93  // Initialise members
94  init_members();
95 
96  // Return
97  return;
98 }
99 
100 
101 /***********************************************************************//**
102  * @brief Type constructor
103  *
104  * @param[in] type Model type.
105  *
106  * Construct an empty sky model of the specified @p type. This constructor
107  * does basically the same than the void constructor with the only difference
108  * that the m_type member of the class is set to the specified @p type
109  * string.
110  ***************************************************************************/
111 GModelSky::GModelSky(const std::string& type) : GModel()
112 {
113  // Initialise members
114  init_members();
115 
116  // Set model type
117  m_type = type;
118 
119  // Return
120  return;
121 }
122 
123 
124 /***********************************************************************//**
125  * @brief XML constructor
126  *
127  * @param[in] xml XML element.
128  *
129  * Constructs a sky model from the model information that is found in the
130  * @p xml element. Only the spatial and spectral information is handled by
131  * this constructor. For the temporal component, a constant model of type
132  * GModelTemporalConst will be allocated.
133  ***************************************************************************/
135 {
136  // Initialise members
137  init_members();
138 
139  // Get pointers on spectrum and spatial model
140  const GXmlElement* spec = xml.element("spectrum", 0);
141  const GXmlElement* spat = xml.element("spatialModel", 0);
142 
143  // Allocate constant
145 
146  // Clone spatial and spectral models
147  m_spatial = xml_spatial(*spat);
148  m_spectral = xml_spectral(*spec);
149  m_temporal = temporal.clone();
150 
151  // Set parameter pointers
152  set_pointers();
153 
154  // Set model type dependent on spatial model type
155  set_type();
156 
157  // Return
158  return;
159 }
160 
161 
162 /***********************************************************************//**
163  * @brief Construct sky model from spatial and spectral XML elements
164  *
165  * @param[in] spatial Spatial XML element.
166  * @param[in] spectral Spectral XML element.
167  *
168  * Constructs a sky model from the @p spatial and @p spectral information
169  * that is found in the respective XML elements. For the temporal component,
170  * a constant model of type GModelTemporalConst will be allocated.
171  ***************************************************************************/
172 GModelSky::GModelSky(const GXmlElement& spatial, const GXmlElement& spectral) :
173  GModel()
174 {
175  // Initialise private members for clean destruction
176  init_members();
177 
178  // Allocate constant
180 
181  // Clone spatial and spectral models
182  m_spatial = xml_spatial(spatial);
183  m_spectral = xml_spectral(spectral);
184  m_temporal = temporal.clone();
185 
186  // Set parameter pointers
187  set_pointers();
188 
189  // Set model type dependent on spatial model type
190  set_type();
191 
192  // Return
193  return;
194 }
195 
196 
197 /***********************************************************************//**
198  * @brief Construct sky model from spatial, spectral and temporal XML elements
199  *
200  * @param[in] spatial Spatial XML element.
201  * @param[in] spectral Spectral XML element.
202  * @param[in] temporal Temporal XML element.
203  *
204  * Constructs a sky model from the @p spatial, @p spectral and @p temporal
205  * information that is found in the respective XML elements.
206  ***************************************************************************/
208  const GXmlElement& spectral,
209  const GXmlElement& temporal) :
210  GModel()
211 {
212  // Initialise private members for clean destruction
213  init_members();
214 
215  // Clone spatial and spectral models
216  m_spatial = xml_spatial(spatial);
217  m_spectral = xml_spectral(spectral);
218  m_temporal = xml_temporal(temporal);
219 
220  // Set parameter pointers
221  set_pointers();
222 
223  // Set model type dependent on spatial model type
224  set_type();
225 
226  // Return
227  return;
228 }
229 
230 
231 /***********************************************************************//**
232  * @brief Construct sky model from spatial and spectral components
233  *
234  * @param[in] spatial Spatial model component.
235  * @param[in] spectral Spectral model component.
236  *
237  * Constructs a sky model from @p spatial and @p spectral model components.
238  * For the temporal component, a constant model of type GModelTemporalConst
239  * will be allocated.
240  ***************************************************************************/
242  const GModelSpectral& spectral) :
243  GModel()
244 {
245  // Initialise members
246  init_members();
247 
248  // Allocate temporal constant model
250 
251  // Clone model components
252  m_spatial = spatial.clone();
253  m_spectral = spectral.clone();
254  m_temporal = temporal.clone();
255 
256  // Set parameter pointers
257  set_pointers();
258 
259  // Set model type dependent on spatial model type
260  set_type();
261 
262  // Return
263  return;
264 }
265 
266 
267 /***********************************************************************//**
268  * @brief Construct sky model from spatial, spectral and temporal components
269  *
270  * @param[in] spatial Spatial model component.
271  * @param[in] spectral Spectral model component.
272  * @param[in] temporal Temporal model component.
273  *
274  * Constructs a sky model from @p spatial, @p spectral and @p temporal model
275  * components.
276  ***************************************************************************/
278  const GModelSpectral& spectral,
279  const GModelTemporal& temporal) :
280  GModel()
281 {
282  // Initialise members
283  init_members();
284 
285  // Clone model components
286  m_spatial = spatial.clone();
287  m_spectral = spectral.clone();
288  m_temporal = temporal.clone();
289 
290  // Set parameter pointers
291  set_pointers();
292 
293  // Set model type dependent on spatial model type
294  set_type();
295 
296  // Return
297  return;
298 }
299 
300 
301 /***********************************************************************//**
302  * @brief Copy constructor
303  *
304  * @param[in] model Sky model.
305  ***************************************************************************/
306 GModelSky::GModelSky(const GModelSky& model) : GModel(model)
307 {
308  // Initialise private members for clean destruction
309  init_members();
310 
311  // Copy members
312  copy_members(model);
313 
314  // Return
315  return;
316 }
317 
318 
319 /***********************************************************************//**
320  * @brief Destructor
321  ***************************************************************************/
323 {
324  // Free members
325  free_members();
326 
327  // Return
328  return;
329 }
330 
331 
332 /*==========================================================================
333  = =
334  = Operators =
335  = =
336  ==========================================================================*/
337 
338 /***********************************************************************//**
339  * @brief Assignment operator
340  *
341  * @param[in] model Sky model.
342  * @return Sky model.
343  ***************************************************************************/
345 {
346  // Execute only if object is not identical
347  if (this != &model) {
348 
349  // Copy base class members
350  this->GModel::operator=(model);
351 
352  // Free members
353  free_members();
354 
355  // Initialise private members for clean destruction
356  init_members();
357 
358  // Copy members
359  copy_members(model);
360 
361  } // endif: object was not identical
362 
363  // Return
364  return *this;
365 }
366 
367 
368 /*==========================================================================
369  = =
370  = Public methods =
371  = =
372  ==========================================================================*/
373 
374 /***********************************************************************//**
375  * @brief Clear sky model
376  *
377  * This method properly resets the sky model to an initial state.
378  ***************************************************************************/
380 {
381  // Free class members
382  free_members();
383  this->GModel::free_members();
384 
385  // Initialise members
386  this->GModel::init_members();
387  init_members();
388 
389  // Return
390  return;
391 }
392 
393 
394 /***********************************************************************//**
395  * @brief Clone sky model
396  *
397  * @return Pointer to deep copy of sky model.
398  ***************************************************************************/
400 {
401  // Clone sky model
402  return new GModelSky(*this);
403 }
404 
405 
406 /***********************************************************************//**
407  * @brief Set spatial model component
408  *
409  * @param[in] spatial Pointer to spatial model component.
410  *
411  * Sets the spatial model component of the model.
412  ***************************************************************************/
413 void GModelSky::spatial(const GModelSpatial* spatial)
414 {
415  // Free spatial model component only if it differs from current
416  // component. This prevents unintential deallocation of the argument
417  if ((m_spatial != NULL) && (m_spatial != spatial)) {
418  delete m_spatial;
419  }
420 
421  // Clone spatial model component if it exists, otherwise set pointer
422  // to NULL
423  m_spatial = (spatial != NULL) ? spatial->clone() : NULL;
424 
425  // Set parameter pointers
426  set_pointers();
427 
428  // Set model type dependent on spatial model type
429  set_type();
430 
431  // Return
432  return;
433 }
434 
435 
436 /***********************************************************************//**
437  * @brief Set spectral model component
438  *
439  * @param[in] spectral Pointer to spectral model component.
440  *
441  * Sets the spectral model component of the model.
442  ***************************************************************************/
443 void GModelSky::spectral(const GModelSpectral* spectral)
444 {
445  // Free spectral model component only if it differs from current
446  // component. This prevents unintential deallocation of the argument
447  if ((m_spectral != NULL) && (m_spectral != spectral)) {
448  delete m_spectral;
449  }
450 
451  // Clone spectral model component if it exists, otherwise set pointer
452  // to NULL
453  m_spectral = (spectral != NULL) ? spectral->clone() : NULL;
454 
455  // Set parameter pointers
456  set_pointers();
457 
458  // Return
459  return;
460 }
461 
462 
463 /***********************************************************************//**
464  * @brief Set temporal model component
465  *
466  * @param[in] temporal Pointer to temporal model component.
467  *
468  * Sets the temporal model component of the model.
469  ***************************************************************************/
470 void GModelSky::temporal(const GModelTemporal* temporal)
471 {
472  // Free temporal model component only if it differs from current
473  // component. This prevents unintential deallocation of the argument
474  if ((m_temporal != NULL) && (m_temporal != temporal)) {
475  delete m_temporal;
476  }
477 
478  // Clone temporal model component if it exists, otherwise set pointer
479  // to NULL
480  m_temporal = (temporal != NULL) ? temporal->clone() : NULL;
481 
482  // Set parameter pointers
483  set_pointers();
484 
485  // Return
486  return;
487 }
488 
489 
490 /***********************************************************************//**
491  * @brief Return value of sky model for a given photon
492  *
493  * @param[in] photon Photon.
494  * @return Value of sky model.
495  *
496  * Returns the value of the sky model for a given @p photon. If no model
497  * components exist the model will return a value of 1.
498  *
499  * @todo We probably should return a value of 0 is no model components exist.
500  * But we may also make sure that the model never has NULL pointers, which
501  * would avoid all the pointer checks.
502  ***************************************************************************/
503 double GModelSky::value(const GPhoton& photon)
504 {
505  // Initialise source model
506  double source = 1.0;
507 
508  // Evaluate source model
509  if (m_spatial != NULL) source *= m_spatial->eval(photon);
510  if (m_spectral != NULL) source *= m_spectral->eval(photon.energy(),
511  photon.time());
512  if (m_temporal != NULL) source *= m_temporal->eval(photon.time());
513 
514  // Return
515  return source;
516 }
517 
518 
519 /***********************************************************************//**
520  * @brief Return parameter gradients of sky model for a given photon
521  *
522  * @param[in] photon Photon.
523  * @return Vector of parameter gradients
524  *
525  * Returns a vector of parameter gradients for the sky model for a given
526  * @p photon. If there are no parameters in the sky model, an empty vector
527  * will be returned.
528  ***************************************************************************/
530 {
531  // Evaluate source model gradients
532  if (m_spatial != NULL) m_spatial->eval(photon, true);
533  if (m_spectral != NULL) m_spectral->eval(photon.energy(), photon.time(),
534  true);
535  if (m_temporal != NULL) m_temporal->eval(photon.time(), true);
536 
537  // Set vector of gradients
539  if (size() > 0) {
540  gradients = GVector(size());
541  for (int i = 0; i < size(); ++i) {
542  gradients[i] = m_pars[i]->factor_gradient();
543  }
544  }
545 
546  // Return gradients
547  return gradients;
548 }
549 
550 
551 /***********************************************************************//**
552  * @brief Evaluate sky model for a given event of an observation
553  *
554  * @param[in] event Observed event.
555  * @param[in] obs Observation.
556  * @param[in] gradients Compute gradients?
557  * @return Value of sky model.
558  *
559  * Evalutes the value of the sky model for an @p event of a specific
560  * observation @p obs.
561  *
562  * If the @p gradients flag is true the method will also compute the
563  * parameter gradients for all model parameters.
564  ***************************************************************************/
565 double GModelSky::eval(const GEvent& event,
566  const GObservation& obs,
567  const bool& gradients) const
568 {
569  // Evaluate model
570  double value = obs.response()->convolve(*this, event, obs, gradients);
571 
572  // If gradients were requested then signal all computed analytical
573  // gradients
574  if (gradients) {
576  }
577 
578  // Return
579  return value;
580 }
581 
582 
583 /***********************************************************************//**
584  * @brief Evaluate sky model for all events of an observation
585  *
586  * @param[in] obs Observation.
587  * @param[out] gradients Pointer to sparse matrix holding model gradients.
588  * @return Values of sky model.
589  *
590  * Evalutes the value of the sky model for an @p event of a specific
591  * observation @p obs.
592  *
593  * If the @p gradients flag is true the method will also compute the
594  * parameter gradients for all model parameters.
595  ***************************************************************************/
597  GMatrixSparse* gradients) const
598 {
599  // Evaluate model
600  GVector values = obs.response()->convolve(*this, obs, gradients);
601 
602  // If gradients were requested then signal all computed analytical
603  // gradients
604  if (gradients) {
606  }
607 
608  // Return
609  return values;
610 }
611 
612 
613 /***********************************************************************//**
614  * @brief Return spatially integrated sky model
615  *
616  * @param[in] obsEng Measured photon energy.
617  * @param[in] obsTime Measured photon arrival time.
618  * @param[in] obs Observation.
619  * @return Spatially integrated sky model.
620  *
621  * Computes
622  *
623  * \f[
624  * N_{\rm pred}(E',t') = \int_{\rm ROI}
625  * P(p',E',t') \, dp' \, dE'
626  * \f]
627  *
628  * of the event probability
629  *
630  * \f[
631  * P(p',E',t') = \int \int \int
632  * S(p,E,t) \times R(p',E',t'|p,E,t) \, dp \, dE \, dt
633  * \f]
634  *
635  * where
636  * \f$S(p,E,t)\f$ is the source model,
637  * \f$R(p',E',t'|p,E,t)\f$ is the instrument response function,
638  * \f$p'\f$ is the measured photon direction,
639  * \f$E'\f$ is the measured photon energy,
640  * \f$t'\f$ is the measured photon arrival time,
641  * \f$p\f$ is the true photon arrival direction,
642  * \f$E\f$ is the true photon energy, and
643  * \f$t\f$ is the true photon arrival time.
644  *
645  * The method calls the GResponse::nroi() method that does the relevant
646  * integration.
647  ***************************************************************************/
648 double GModelSky::npred(const GEnergy& obsEng,
649  const GTime& obsTime,
650  const GObservation& obs) const
651 {
652  // Initialise result
653  double npred = 0.0;
654 
655  // Continue only if model is valid)
656  if (valid_model()) {
657 
658  // Compute Nroi
659  npred = obs.response()->nroi(*this, obsEng, obsTime, obs);
660 
661  // Compile option: Check for NaN/Inf
662  #if defined(G_NAN_CHECK)
663  if (gammalib::is_notanumber(npred) || gammalib::is_infinite(npred)) {
664  std::cout << "*** ERROR: GModelSky::npred:";
665  std::cout << " NaN/Inf encountered";
666  std::cout << " (npred=" << npred;
667  std::cout << ", obsEng=" << obsEng;
668  std::cout << ", obsTime=" << obsTime;
669  std::cout << ")" << std::endl;
670  }
671  #endif
672 
673  } // endif: model was valid
674 
675  // Return npred
676  return npred;
677 }
678 
679 
680 /***********************************************************************//**
681  * @brief Read sky model from XML element
682  *
683  * @param[in] xml XML element.
684  *
685  * Reads the sky model from an XML element. The XML element is expected to
686  * respect the following format:
687  *
688  * <source name=".." type=".." instrument=".." id=".." ts="..">
689  * <spectrum type="..">
690  * ..
691  * </spectrum>
692  * <spatialModel type="..">
693  * ..
694  * </spatialModel>
695  * <temporal type="..">
696  * ..
697  * </temporal>
698  * </source>
699  *
700  * The temporal element is optional. In no temporal element is specified a
701  * constant component with unity normalization will be assumed.
702  *
703  * Optionally, the model may also contain scale parameters following the
704  * format:
705  *
706  * <source name=".." type=".." instrument=".." id=".." ts="..">
707  * <spectrum type="..">
708  * ..
709  * </spectrum>
710  * <spatialModel type="..">
711  * ..
712  * </spatialModel>
713  * <temporal type="..">
714  * ..
715  * </temporal>
716  * <scaling>
717  * <instrument name=".." scale="1.0" min="0.1" max="10.0" value="1.0" free="0"/>
718  * <instrument name=".." scale="1.0" min="0.1" max="10.0" value="0.5" free="0"/>
719  * </scaling>
720  * </source>
721  *
722  * (see GModel::read_scales() for more information on instrument scales).
723  ***************************************************************************/
724 void GModelSky::read(const GXmlElement& xml)
725 {
726  // Clear sky model
727  clear();
728 
729  // Get pointers on spectrum and spatial model
730  const GXmlElement* spec = xml.element("spectrum", 0);
731  const GXmlElement* spat = xml.element("spatialModel", 0);
732 
733  // Set spatial and spectral models
734  m_spatial = xml_spatial(*spat);
735  m_spectral = xml_spectral(*spec);
736 
737  // Handle optional temporal model
738  if (xml.elements("temporal") > 0) {
739  const GXmlElement* temp = xml.element("temporal", 0);
740  m_temporal = xml_temporal(*temp);
741  }
742  else {
744  m_temporal = temporal.clone();
745  }
746 
747  // Read model attributes
748  read_attributes(xml);
749 
750  // Set parameter pointers
751  set_pointers();
752 
753  // Set model type dependent on spatial model type
754  set_type();
755 
756  // Return
757  return;
758 }
759 
760 
761 /***********************************************************************//**
762  * @brief Write model into XML element
763  *
764  * @param[in] xml Source library.
765  *
766  * Writes the sky model into an XML source library. The format written to
767  * the @p xml element is as follows:
768  *
769  * <source name=".." type=".." instrument=".." id=".." ts="..">
770  * <spectrum type="..">
771  * ..
772  * </spectrum>
773  * <spatialModel type="..">
774  * ..
775  * </spatialModel>
776  * <temporal type="..">
777  * ..
778  * </temporal>
779  * <scaling>
780  * <instrument name=".." scale="1.0" min="0.1" max="10.0" value="1.0" free="0"/>
781  * <instrument name=".." scale="1.0" min="0.1" max="10.0" value="0.5" free="0"/>
782  * </scaling>
783  * </source>
784  *
785  * For compatibility reasons the temporal element will only be written if it
786  * is a non-constant component or a constant component with a normalization
787  * that differs from unity.
788  *
789  * The scaling element will only be written optionally in case that instrument
790  * dependent scaling factors exist (see GModel::write_scales() for more
791  * information on instrument scales).
792  ***************************************************************************/
794 {
795  // Initialise pointer on source
796  GXmlElement* src = NULL;
797 
798  // Search corresponding source
799  int n = xml.elements("source");
800  for (int k = 0; k < n; ++k) {
801  GXmlElement* element = static_cast<GXmlElement*>(xml.element("source", k));
802  if (element->attribute("name") == name()) {
803  src = element;
804  break;
805  }
806  }
807 
808  // If the temporal model is not a constant with unit normalization then
809  // set cons to a NULL pointer
810  GModelTemporalConst* cons = dynamic_cast<GModelTemporalConst*>(temporal());
811  if (cons != NULL) {
812  if (cons->norm() != 1.0) {
813  cons = NULL;
814  }
815  }
816 
817  // If no source with corresponding name was found then append one
818  if (src == NULL) {
819  src = xml.append("source");
820  if (spectral() != NULL) src->append(GXmlElement("spectrum"));
821  if (spatial() != NULL) src->append(GXmlElement("spatialModel"));
822  if (temporal() != NULL && cons == NULL) src->append(GXmlElement("temporal"));
823  }
824 
825  // Write spectral model
826  if (spectral() != NULL) {
827  GXmlElement* spec = src->element("spectrum", 0);
828  spectral()->write(*spec);
829  }
830 
831  // Write spatial model
832  if (spatial() != NULL) {
833  GXmlElement* spat = src->element("spatialModel", 0);
834  spatial()->write(*spat);
835  }
836 
837  // Write temporal model (only if not a constant with unit normalization
838  // factor)
839  if (temporal() != NULL && cons == NULL) {
840  GXmlElement* temp = src->element("temporal", 0);
841  temporal()->write(*temp);
842  }
843 
844  // Write model attributes
845  write_attributes(*src);
846 
847  // Return
848  return;
849 }
850 
851 
852 /***********************************************************************//**
853  * @brief Return simulated list of photons
854  *
855  * @param[in] area Simulation surface area (cm2).
856  * @param[in] dir Centre of simulation cone.
857  * @param[in] radius Radius of simulation cone (deg).
858  * @param[in] emin Minimum photon energy.
859  * @param[in] emax Maximum photon energy.
860  * @param[in] tmin Minimum photon arrival time.
861  * @param[in] tmax Maximum photon arrival time.
862  * @param[in,out] ran Random number generator.
863  * @return List of photons
864  *
865  * Returns a list of photons that has been derived by Monte Carlo simulation
866  * from the model. A simulation region is define by specification of
867  * - a simulation cone, which is a circular region on the sky defined by
868  * a centre direction @p dir and a @p radius,
869  * - an energy range [@p emin, @p emax], and
870  * - a time interval [@p tmin, @p tmax].
871  *
872  * Only photons with parameters in the simulation region will be returned
873  * by the method.
874  *
875  * The simulation cone may eventually cover the entire sky (by setting
876  * the radius to 180 degrees), yet simulations will be more efficient if
877  * only the sky region will be simulated that is actually observed by the
878  * telescope.
879  *
880  * @todo Implement unique model ID to assign as Monte Carlo ID
881  ***************************************************************************/
882 GPhotons GModelSky::mc(const double& area,
883  const GSkyDir& dir, const double& radius,
884  const GEnergy& emin, const GEnergy& emax,
885  const GTime& tmin, const GTime& tmax,
886  GRan& ran) const
887 {
888  // Allocate photons
889  GPhotons photons;
890 
891  // Continue only if model is valid)
892  if (valid_model()) {
893 
894  // Determine the spatial model normalization within the simulation
895  // cone and check whether the model will produce any photons in that
896  // cone.
897  double norm = m_spatial->mc_norm(dir, radius);
898  bool use_model = (norm > 0.0) ? true : false;
899 
900  // Continue only if model overlaps with simulation region
901  if (use_model) {
902 
903  // Initialise de-allocation flag
904  bool free_spectral = false;
905 
906  // Set pointer to spectral model
908 
909  // If the spectral model is a diffuse cube then create a node
910  // function spectral model that is the product of the diffuse
911  // cube node function and the spectral model evaluated at the
912  // energies of the node function
914  dynamic_cast<GModelSpatialDiffuseCube*>(m_spatial);
915  if (cube != NULL) {
916 
917  // Set MC cone
918  cube->mc_cone(GSkyRegionCircle(dir, radius));
919 
920  // Allocate node function to replace the spectral component
921  GModelSpectralNodes* nodes =
922  new GModelSpectralNodes(cube->spectrum());
923  for (int i = 0; i < nodes->nodes(); ++i) {
924  GEnergy energy = nodes->energy(i);
925  double intensity = nodes->intensity(i);
926  double value = m_spectral->eval(energy, tmin);
927  nodes->intensity(i, value*intensity);
928  }
929 
930  // Signal that node function needs to be de-allocated later
931  free_spectral = true;
932 
933  // Set the spectral model pointer to the node function
934  spectral = nodes;
935 
936  // Kluge: if there are no nodes then the spectral->flux method
937  // will throw an exception. We therefore set here the use_model
938  // flag to false in case that there are no spectral nodes
939  if (nodes->nodes() == 0) {
940  use_model = false;
941  }
942 
943  } // endif: spatial model was a diffuse cube
944 
945  // Compute flux within [emin, emax] in model from spectral
946  // component (units: ph/cm2/s)
947  double flux = (use_model) ? spectral->flux(emin, emax) : 0.0;
948 
949  // Derive expecting counting rate within simulation surface
950  // (units: ph/s)
951  double rate = flux * area * norm;
952 
953  // Debug option: dump rate
954  #if defined(G_DUMP_MC)
955  std::cout << "GModelSky::mc(\"" << name() << "\": ";
956  std::cout << "flux=" << flux << " ph/cm2/s, ";
957  std::cout << "rate=" << rate << " ph/s, ";
958  std::cout << "norm=" << norm << ")" << std::endl;
959  #endif
960 
961  // Continue only if rate is positive
962  if (rate > 0.0) {
963 
964  // Get photon arrival times from temporal model
965  GTimes times = m_temporal->mc(rate, tmin, tmax, ran);
966 
967  // Debug option: dump number of times
968  #if defined(G_DUMP_MC_DETAIL)
969  std::cout << " Times=" << times.size() << std::endl;
970  #endif
971 
972  // Reserve space for photons
973  if (times.size() > 0) {
974  photons.reserve(times.size());
975  }
976 
977  // Loop over photons
978  for (int i = 0; i < times.size(); ++i) {
979 
980  // Debug option: dump photon index
981  #if defined(G_DUMP_MC_DETAIL)
982  std::cout << " Photon=" << i << std::endl;
983  #endif
984 
985  // Allocate photon
986  GPhoton photon;
987 
988  // Set photon arrival time
989  photon.time(times[i]);
990 
991  // Debug option: dump time
992  #if defined(G_DUMP_MC_DETAIL)
993  std::cout << " Time=" << times[i] << std::endl;
994  #endif
995 
996  // Set photon energy. If an invalid_return_value exception
997  // occurs the energy returned by the spectral Monte Carlo
998  // method is invalid and the photon is skipped.
999  try {
1000  photon.energy(spectral->mc(emin, emax, photon.time(),
1001  ran));
1002  }
1004  continue;
1005  }
1006 
1007  // Debug option: dump energy
1008  #if defined(G_DUMP_MC_DETAIL)
1009  std::cout << " Energy=" << photon.energy() << std::endl;
1010  #endif
1011 
1012  // Set incident photon direction. If an invalid_return_value
1013  // exception occurs the sky direction returned by the
1014  // spatial Monte Carlo method is invalid and the photon is
1015  // skipped.
1016  try {
1017  photon.dir(m_spatial->mc(photon.energy(), photon.time(),
1018  ran));
1019  }
1021  continue;
1022  }
1023 
1024  // Debug option: dump direction
1025  #if defined(G_DUMP_MC_DETAIL)
1026  std::cout << " Direction=" << photon.dir() << std::endl;
1027  #endif
1028 
1029  // Append photon
1030  if (dir.dist_deg(photon.dir()) <= radius) {
1031  photons.append(photon);
1032  }
1033 
1034  } // endfor: looped over photons
1035 
1036  } // endif: rate was positive
1037 
1038  // Free spectral model if required
1039  if (free_spectral) delete spectral;
1040 
1041  } // endif: model was used
1042 
1043  } // endif: model was valid
1044 
1045  // Return photon list
1046  return photons;
1047 }
1048 
1049 
1050 /***********************************************************************//**
1051  * @brief Return sky model photon flux
1052  *
1053  * @param[in] emin Minimum energy.
1054  * @param[in] emax Maximum energy.
1055  * @return Sky model photon flux (cm\f$^{-2}\f$ s\f$^{-1}\f$)
1056  *
1057  * Returns the spatially integrate photon flux in the sky model.
1058  *
1059  * @todo Take spatial photon flux normalisation into account
1060  ***************************************************************************/
1061 double GModelSky::flux(const GEnergy& emin, const GEnergy& emax) const
1062 {
1063  // Initialise flux
1064  double flux = 0.0;
1065 
1066  // Get pointer on diffuse cube spatial model component
1067  const GModelSpatialDiffuseCube* cube =
1068  dynamic_cast<const GModelSpatialDiffuseCube*>(spatial());
1069 
1070  // If spatial model is a diffuse cube then compute the flux by
1071  // multiplying the spectrum of the spatial and spectral components
1072  // using a node function
1073  if (cube != NULL) {
1074 
1075  // Get spectral nodes
1076  GModelSpectralNodes nodes = cube->spectrum();
1077 
1078  // Multiply spectral nodes with spectral model
1079  for (int i = 0; i < nodes.nodes(); ++i) {
1080 
1081  // Get energy and intensity
1082  GEnergy energy = nodes.energy(i);
1083  double intensity = nodes.intensity(i);
1084 
1085  // Multiply intensity with spectral model
1086  intensity *= spectral()->eval(energy);
1087 
1088  // Store intensity
1089  nodes.intensity(i, intensity);
1090 
1091  } // endfor: multiplied spectral nodes with spectral model
1092 
1093  // Compute flux
1094  flux = nodes.flux(emin, emax);
1095 
1096  } // endif: spatial model was a diffuse cube
1097 
1098  // ... otherwise assume that the spatial model is unity and compute
1099  // flux from the spectral component
1100  else {
1101  flux = spectral()->flux(emin, emax);
1102  }
1103 
1104  // Return flux
1105  return flux;
1106 }
1107 
1108 
1109 /***********************************************************************//**
1110  * @brief Return sky model photon flux within sky region
1111  *
1112  * @param[in] region Sky region.
1113  * @param[in] emin Minimum energy.
1114  * @param[in] emax Maximum energy.
1115  * @return Sky model photon flux (cm\f$^{-2}\f$ s\f$^{-1}\f$)
1116  *
1117  * Returns the photon flux of the sky model within a given sky region.
1118  ***************************************************************************/
1119 double GModelSky::flux(const GSkyRegion& region,
1120  const GEnergy& emin,
1121  const GEnergy& emax) const
1122 {
1123  // Initialise flux
1124  double flux = 0.0;
1125 
1126  // Signal for which models the spatial component depends on energy
1127  bool dependence = (spatial()->classname() == "GModelSpatialDiffuseCube") ||
1128  (spatial()->classname() == "GModelSpatialComposite");
1129 
1130  // If the spatial component depends on energy then use
1131  if (dependence) {
1132 
1133  // Get boundaries in MeV
1134  double mev_min = emin.MeV();
1135  double mev_max = emax.MeV();
1136 
1137  // Continue only if valid
1138  if (mev_max > mev_min) {
1139 
1140  // Setup integration function
1141  GModelSky::flux_kern integrand(this, &region);
1142  GIntegral integral(&integrand);
1143 
1144  // Do Romberg integration
1145  flux = integral.romberg(std::log(mev_min), std::log(mev_max));
1146 
1147  } // endif: energy interval was valid
1148 
1149  } // endif: spatial model was a diffuse cube
1150 
1151  // ... otherwise the spatial model does not depend on energy, hence
1152  // the flux is computed as the product between the spatial and
1153  // spectral flux. Note that this correctly takes the spatial
1154  // normalisation into account since the flux() method does this.
1155  else {
1156  flux = spatial()->flux(region) * spectral()->flux(emin, emax);
1157  }
1158 
1159  // Return flux
1160  return flux;
1161 }
1162 
1163 
1164 /***********************************************************************//**
1165  * @brief Return sky model energy flux
1166  *
1167  * @param[in] emin Minimum energy.
1168  * @param[in] emax Maximum energy.
1169  * @return Sky model energy flux (erg cm\f$^{-2}\f$ s\f$^{-1}\f$)
1170  *
1171  * Returns the spatially integrate energy flux in the sky model.
1172  ***************************************************************************/
1173 double GModelSky::eflux(const GEnergy& emin, const GEnergy& emax) const
1174 {
1175  // Initialise energy flux
1176  double eflux = 0.0;
1177 
1178  // Get pointer on diffuse cube spatial model component
1179  const GModelSpatialDiffuseCube* cube =
1180  dynamic_cast<const GModelSpatialDiffuseCube*>(spatial());
1181 
1182  // If spatial model is a diffuse cube then compute the flux by
1183  // multiplying the spectrum of the spatial and spectral components
1184  // using a node function
1185  if (cube != NULL) {
1186 
1187  // Get spectral nodes
1188  GModelSpectralNodes nodes = cube->spectrum();
1189 
1190  // Multiply spectral nodes with spectral model
1191  for (int i = 0; i < nodes.nodes(); ++i) {
1192 
1193  // Get energy and intensity
1194  GEnergy energy = nodes.energy(i);
1195  double intensity = nodes.intensity(i);
1196 
1197  // Multiply intensity with spectral model
1198  intensity *= spectral()->eval(energy);
1199 
1200  // Store intensity
1201  nodes.intensity(i, intensity);
1202 
1203  } // endfor: multiplied spectral nodes with spectral model
1204 
1205  // Compute energy flux
1206  eflux = nodes.eflux(emin, emax);
1207 
1208  } // endif: spatial model was a diffuse cube
1209 
1210  // ... otherwise assume that the spatial model is unity and compute
1211  // energy flux from the spectral component
1212  else {
1213  eflux = spectral()->eflux(emin, emax);
1214  }
1215 
1216  // Return energy flux
1217  return eflux;
1218 }
1219 
1220 
1221 /***********************************************************************//**
1222  * @brief Return sky model energy flux within sky region
1223  *
1224  * @param[in] region Sky region.
1225  * @param[in] emin Minimum energy.
1226  * @param[in] emax Maximum energy.
1227  * @return Sky model energy flux (erg cm\f$^{-2}\f$ s\f$^{-1}\f$)
1228  *
1229  * Returns the energy flux of the sky model within a given sky region. If
1230  * the spatial model is energy dependent the method integrates the flux over
1231  * the specified energy range. Otherwise the spatially integrated flux is
1232  * multiplied by the spectral flux.
1233  ***************************************************************************/
1234 double GModelSky::eflux(const GSkyRegion& region,
1235  const GEnergy& emin,
1236  const GEnergy& emax) const
1237 {
1238  // Initialise flux
1239  double eflux = 0.0;
1240 
1241  // Signal for which models the spatial component depends on energy
1242  bool dependence = (spatial()->classname() == "GModelSpatialDiffuseCube") ||
1243  (spatial()->classname() == "GModelSpatialComposite");
1244 
1245  // If the spatial component depends on energy then use
1246  if (dependence) {
1247 
1248  // Get boundaries in MeV
1249  double mev_min = emin.MeV();
1250  double mev_max = emax.MeV();
1251 
1252  // Continue only if valid
1253  if (mev_max > mev_min) {
1254 
1255  // Setup integration function
1256  GModelSky::eflux_kern integrand(this, &region);
1257  GIntegral integral(&integrand);
1258 
1259  // Do Romberg integration
1260  eflux = integral.romberg(std::log(mev_min), std::log(mev_max));
1261 
1262  } // endif: energy interval was valid
1263 
1264  } // endif: spatial model was a diffuse cube
1265 
1266  // ... otherwise the spatial model does not depend on energy, hence
1267  // the flux is computed as the product between the spatial and
1268  // spectral flux. Note that this correctly takes the spatial
1269  // normalisation into account since the flux() method does this.
1270  else {
1271  eflux = spatial()->flux(region) * spectral()->eflux(emin, emax);
1272  }
1273 
1274  // Return flux
1275  return eflux;
1276 }
1277 
1278 
1279 /***********************************************************************//**
1280  * @brief Return sky model photon flux error
1281  *
1282  * @param[in] emin Minimum energy.
1283  * @param[in] emax Maximum energy.
1284  * @return Sky model photon flux error (cm\f$^{-2}\f$ s\f$^{-1}\f$)
1285  *
1286  * Returns the spatially integrate photon flux error in the sky model. The
1287  * photon flux error is computed using Gaussian error propagation
1288  *
1289  * \f[
1290  * \Delta F_{ph} = \sqrt{
1291  * \sum_i \left( \frac{\partial F_{ph}}{\partial p_i} \right)^2
1292  * \left( \Delta p_i \right)^2 }
1293  * \f]
1294  *
1295  * where
1296  * \f$F_{ph}\f$ is the sky model photon flux and
1297  * \f$p_i\f$ are the sky model parameters.
1298  ***************************************************************************/
1299 double GModelSky::flux_error(const GEnergy& emin, const GEnergy& emax) const
1300 {
1301  // Initialise flux error
1302  double flux_error = 0.0;
1303 
1304  // Loop over all model parameters
1305  for (int ipar = 0; ipar < size(); ++ipar) {
1306 
1307  // Get non-const model pointer
1308  GModelPar* par = const_cast<GModelPar*>(&(*this)[ipar]);
1309 
1310  // Consider only free parameters
1311  if (par->is_free()) {
1312 
1313  // Save current model parameter
1314  GModelPar current = *par;
1315 
1316  // Autoscale model parameter
1317  par->autoscale();
1318 
1319  // Get actual parameter value
1320  double x = par->factor_value();
1321 
1322  // Set fixed step size for computation of derivative.
1323  const double step_size = 0.0002;
1324  double h = step_size;
1325 
1326  // Re-adjust the step-size h in case that the initial step size is
1327  // larger than the allowed parameter range
1328  if (par->has_min() && par->has_max()) {
1329  double par_h = par->factor_max() - par->factor_min();
1330  if (par_h < h) {
1331  h = par_h;
1332  }
1333  }
1334 
1335  // Initialise derivative
1336  double derivative = 0.0;
1337 
1338  // Continue only if step size is positive
1339  if (h > 0.0) {
1340 
1341  // If we are too close to the minimum boundary use a right
1342  // sided difference ...
1343  if (par->has_min() && ((x-par->factor_min()) < h)) {
1344 
1345  // Compute fs1
1346  par->factor_value(x);
1347  double fs1 = flux(emin, emax);
1348 
1349  // Compute fs2
1350  par->factor_value(x+h);
1351  double fs2 = flux(emin, emax);
1352 
1353  // Compute derivative
1354  derivative = (fs1 - fs2) / h;
1355 
1356  }
1357 
1358  // ... otherwise if we are too close to the maximum boundary use
1359  // a left sided difference ...
1360  else if (par->has_max() && ((par->factor_max()-x) < h)) {
1361 
1362  // Compute fs1
1363  par->factor_value(x);
1364  double fs1 = flux(emin, emax);
1365 
1366  // Compute fs2
1367  par->factor_value(x-h);
1368  double fs2 = flux(emin, emax);
1369 
1370  // Compute derivative
1371  derivative = (fs1 - fs2) / h;
1372 
1373  }
1374 
1375  // ... otherwise use a symmetric difference
1376  else {
1377 
1378  // Compute fs1
1379  par->factor_value(x+h);
1380  double fs1 = flux(emin, emax);
1381 
1382  // Compute fs2
1383  par->factor_value(x-h);
1384  double fs2 = flux(emin, emax);
1385 
1386  // Compute derivative
1387  derivative = 0.5 * (fs1 - fs2) / h;
1388 
1389  }
1390 
1391  } // endif: step size was positive
1392 
1393  // Add contribution
1394  if (derivative > 0.0) {
1395  double factor = derivative * par->factor_error();
1396  flux_error += factor * factor;
1397  }
1398 
1399  // Restore current model parameter
1400  *par = current;
1401 
1402  } // endif: parameter was free
1403 
1404  } // endfor: looped over all parameters
1405 
1406  // Take square root if required
1407  if (flux_error > 0.0) {
1408  flux_error = std::sqrt(flux_error);
1409  }
1410 
1411  // Return flux error
1412  return flux_error;
1413 }
1414 
1415 
1416 /***********************************************************************//**
1417  * @brief Return sky model photon flux error within sky region
1418  *
1419  * @param[in] region Sky region.
1420  * @param[in] emin Minimum energy.
1421  * @param[in] emax Maximum energy.
1422  * @return Sky model photon flux error (cm\f$^{-2}\f$ s\f$^{-1}\f$)
1423  *
1424  * Returns the spatially integrate photon flux error in the sky model. The
1425  * photon flux error is computed using Gaussian error propagation
1426  *
1427  * \f[
1428  * \Delta F_{ph} = \sqrt{
1429  * \sum_i \left( \frac{\partial F_{ph}}{\partial p_i} \right)^2
1430  * \left( \Delta p_i \right)^2 }
1431  * \f]
1432  *
1433  * where
1434  * \f$F_{ph}\f$ is the sky model photon flux and
1435  * \f$p_i\f$ are the sky model parameters.
1436  ***************************************************************************/
1437 double GModelSky::flux_error(const GSkyRegion& region,
1438  const GEnergy& emin,
1439  const GEnergy& emax) const
1440 {
1441  // Initialise flux error
1442  double flux_error = 0.0;
1443 
1444  // Loop over all model parameters
1445  for (int ipar = 0; ipar < size(); ++ipar) {
1446 
1447  // Get non-const model pointer
1448  GModelPar* par = const_cast<GModelPar*>(&(*this)[ipar]);
1449 
1450  // Consider only free parameters
1451  if (par->is_free()) {
1452 
1453  // Save current model parameter
1454  GModelPar current = *par;
1455 
1456  // Autoscale model parameter
1457  par->autoscale();
1458 
1459  // Get actual parameter value
1460  double x = par->factor_value();
1461 
1462  // Set fixed step size for computation of derivative.
1463  const double step_size = 0.0002;
1464  double h = step_size;
1465 
1466  // Re-adjust the step-size h in case that the initial step size is
1467  // larger than the allowed parameter range
1468  if (par->has_min() && par->has_max()) {
1469  double par_h = par->factor_max() - par->factor_min();
1470  if (par_h < h) {
1471  h = par_h;
1472  }
1473  }
1474 
1475  // Initialise derivative
1476  double derivative = 0.0;
1477 
1478  // Continue only if step size is positive
1479  if (h > 0.0) {
1480 
1481  // If we are too close to the minimum boundary use a right
1482  // sided difference ...
1483  if (par->has_min() && ((x-par->factor_min()) < h)) {
1484 
1485  // Compute fs1
1486  par->factor_value(x);
1487  double fs1 = flux(region, emin, emax);
1488 
1489  // Compute fs2
1490  par->factor_value(x+h);
1491  double fs2 = flux(region, emin, emax);
1492 
1493  // Compute derivative
1494  derivative = (fs1 - fs2) / h;
1495 
1496  }
1497 
1498  // ... otherwise if we are too close to the maximum boundary use
1499  // a left sided difference ...
1500  else if (par->has_max() && ((par->factor_max()-x) < h)) {
1501 
1502  // Compute fs1
1503  par->factor_value(x);
1504  double fs1 = flux(region, emin, emax);
1505 
1506  // Compute fs2
1507  par->factor_value(x-h);
1508  double fs2 = flux(region, emin, emax);
1509 
1510  // Compute derivative
1511  derivative = (fs1 - fs2) / h;
1512 
1513  }
1514 
1515  // ... otherwise use a symmetric difference
1516  else {
1517 
1518  // Compute fs1
1519  par->factor_value(x+h);
1520  double fs1 = flux(region, emin, emax);
1521 
1522  // Compute fs2
1523  par->factor_value(x-h);
1524  double fs2 = flux(region, emin, emax);
1525 
1526  // Compute derivative
1527  derivative = 0.5 * (fs1 - fs2) / h;
1528 
1529  }
1530 
1531  } // endif: step size was positive
1532 
1533  // Add contribution
1534  if (derivative > 0.0) {
1535  double factor = derivative * par->factor_error();
1536  flux_error += factor * factor;
1537  }
1538 
1539  // Restore current model parameter
1540  *par = current;
1541 
1542  } // endif: parameter was free
1543 
1544  } // endfor: looped over all parameters
1545 
1546  // Take square root if required
1547  if (flux_error > 0.0) {
1548  flux_error = std::sqrt(flux_error);
1549  }
1550  // Return flux error
1551  return flux_error;
1552 }
1553 
1554 
1555 /***********************************************************************//**
1556  * @brief Return sky model energy flux error
1557  *
1558  * @param[in] emin Minimum energy.
1559  * @param[in] emax Maximum energy.
1560  * @return Sky model energy flux error (erg cm\f$^{-2}\f$ s\f$^{-1}\f$)
1561  *
1562  * Returns the spatially integrate energy flux error in the sky model. The
1563  * energy flux error is computed using Gaussian error propagation
1564  *
1565  * \f[
1566  * \Delta F_E = \sqrt{
1567  * \sum_i \left( \frac{\partial F_E}{\partial p_i} \right)^2
1568  * \left( \Delta p_i \right)^2 }
1569  * \f]
1570  *
1571  * where
1572  * \f$F_E\f$ is the sky model energy flux and
1573  * \f$p_i\f$ are the sky model parameters.
1574  ***************************************************************************/
1575 double GModelSky::eflux_error(const GEnergy& emin, const GEnergy& emax) const
1576 {
1577  // Initialise flux error
1578  double eflux_error = 0.0;
1579 
1580  // Loop over all model parameters
1581  for (int ipar = 0; ipar < size(); ++ipar) {
1582 
1583  // Get non-const model pointer
1584  GModelPar* par = const_cast<GModelPar*>(&(*this)[ipar]);
1585 
1586  // Consider only free parameters
1587  if (par->is_free()) {
1588 
1589  // Save current model parameter
1590  GModelPar current = *par;
1591 
1592  // Autoscale model parameter
1593  par->autoscale();
1594 
1595  // Get actual parameter value
1596  double x = par->factor_value();
1597 
1598  // Set fixed step size for computation of derivative.
1599  const double step_size = 0.0002;
1600  double h = step_size;
1601 
1602  // Re-adjust the step-size h in case that the initial step size is
1603  // larger than the allowed parameter range
1604  if (par->has_min() && par->has_max()) {
1605  double par_h = par->factor_max() - par->factor_min();
1606  if (par_h < h) {
1607  h = par_h;
1608  }
1609  }
1610 
1611  // Initialise derivative
1612  double derivative = 0.0;
1613 
1614  // Continue only if step size is positive
1615  if (h > 0.0) {
1616 
1617  // If we are too close to the minimum boundary use a right
1618  // sided difference ...
1619  if (par->has_min() && ((x-par->factor_min()) < h)) {
1620 
1621  // Compute fs1
1622  par->factor_value(x);
1623  double fs1 = eflux(emin, emax);
1624 
1625  // Compute fs2
1626  par->factor_value(x+h);
1627  double fs2 = eflux(emin, emax);
1628 
1629  // Compute derivative
1630  derivative = (fs1 - fs2) / h;
1631 
1632  }
1633 
1634  // ... otherwise if we are too close to the maximum boundary use
1635  // a left sided difference ...
1636  else if (par->has_max() && ((par->factor_max()-x) < h)) {
1637 
1638  // Compute fs1
1639  par->factor_value(x);
1640  double fs1 = eflux(emin, emax);
1641 
1642  // Compute fs2
1643  par->factor_value(x-h);
1644  double fs2 = eflux(emin, emax);
1645 
1646  // Compute derivative
1647  derivative = (fs1 - fs2) / h;
1648 
1649  }
1650 
1651  // ... otherwise use a symmetric difference
1652  else {
1653 
1654  // Compute fs1
1655  par->factor_value(x+h);
1656  double fs1 = eflux(emin, emax);
1657 
1658  // Compute fs2
1659  par->factor_value(x-h);
1660  double fs2 = eflux(emin, emax);
1661 
1662  // Compute derivative
1663  derivative = 0.5 * (fs1 - fs2) / h;
1664 
1665  }
1666 
1667  } // endif: step size was positive
1668 
1669  // Add contribution
1670  if (derivative > 0.0) {
1671  double factor = derivative * par->factor_error();
1672  eflux_error += factor * factor;
1673  }
1674 
1675  // Restore current model parameter
1676  *par = current;
1677 
1678  } // endif: parameter was free
1679 
1680  } // endfor: looped over all parameters
1681 
1682  // Take square root if required
1683  if (eflux_error > 0.0) {
1684  eflux_error = std::sqrt(eflux_error);
1685  }
1686 
1687  // Return flux error
1688  return eflux_error;
1689 }
1690 
1691 
1692 /***********************************************************************//**
1693  * @brief Return sky model energy flux error within sky region
1694  *
1695  * @param[in] region Sky region.
1696  * @param[in] emin Minimum energy.
1697  * @param[in] emax Maximum energy.
1698  * @return Sky model energy flux error (erg cm\f$^{-2}\f$ s\f$^{-1}\f$)
1699  *
1700  * Returns the spatially integrate energy flux error in the sky model. The
1701  * energy flux error is computed using Gaussian error propagation
1702  *
1703  * \f[
1704  * \Delta F_E = \sqrt{
1705  * \sum_i \left( \frac{\partial F_E}{\partial p_i} \right)^2
1706  * \left( \Delta p_i \right)^2 }
1707  * \f]
1708  *
1709  * where
1710  * \f$F_E\f$ is the sky model energy flux and
1711  * \f$p_i\f$ are the sky model parameters.
1712  ***************************************************************************/
1713 double GModelSky::eflux_error(const GSkyRegion& region,
1714  const GEnergy& emin,
1715  const GEnergy& emax) const
1716 {
1717  // Initialise flux error
1718  double eflux_error = 0.0;
1719 
1720  // Loop over all model parameters
1721  for (int ipar = 0; ipar < size(); ++ipar) {
1722 
1723  // Get non-const model pointer
1724  GModelPar* par = const_cast<GModelPar*>(&(*this)[ipar]);
1725 
1726  // Consider only free parameters
1727  if (par->is_free()) {
1728 
1729  // Save current model parameter
1730  GModelPar current = *par;
1731 
1732  // Autoscale model parameter
1733  par->autoscale();
1734 
1735  // Get actual parameter value
1736  double x = par->factor_value();
1737 
1738  // Set fixed step size for computation of derivative.
1739  const double step_size = 0.0002;
1740  double h = step_size;
1741 
1742  // Re-adjust the step-size h in case that the initial step size is
1743  // larger than the allowed parameter range
1744  if (par->has_min() && par->has_max()) {
1745  double par_h = par->factor_max() - par->factor_min();
1746  if (par_h < h) {
1747  h = par_h;
1748  }
1749  }
1750 
1751  // Initialise derivative
1752  double derivative = 0.0;
1753 
1754  // Continue only if step size is positive
1755  if (h > 0.0) {
1756 
1757  // If we are too close to the minimum boundary use a right
1758  // sided difference ...
1759  if (par->has_min() && ((x-par->factor_min()) < h)) {
1760 
1761  // Compute fs1
1762  par->factor_value(x);
1763  double fs1 = eflux(region, emin, emax);
1764 
1765  // Compute fs2
1766  par->factor_value(x+h);
1767  double fs2 = eflux(region, emin, emax);
1768 
1769  // Compute derivative
1770  derivative = (fs1 - fs2) / h;
1771 
1772  }
1773 
1774  // ... otherwise if we are too close to the maximum boundary use
1775  // a left sided difference ...
1776  else if (par->has_max() && ((par->factor_max()-x) < h)) {
1777 
1778  // Compute fs1
1779  par->factor_value(x);
1780  double fs1 = eflux(region, emin, emax);
1781 
1782  // Compute fs2
1783  par->factor_value(x-h);
1784  double fs2 = eflux(region, emin, emax);
1785 
1786  // Compute derivative
1787  derivative = (fs1 - fs2) / h;
1788 
1789  }
1790 
1791  // ... otherwise use a symmetric difference
1792  else {
1793 
1794  // Compute fs1
1795  par->factor_value(x+h);
1796  double fs1 = eflux(region, emin, emax);
1797 
1798  // Compute fs2
1799  par->factor_value(x-h);
1800  double fs2 = eflux(region, emin, emax);
1801 
1802  // Compute derivative
1803  derivative = 0.5 * (fs1 - fs2) / h;
1804 
1805  }
1806 
1807  } // endif: step size was positive
1808 
1809  // Add contribution
1810  if (derivative > 0.0) {
1811  double factor = derivative * par->factor_error();
1812  eflux_error += factor * factor;
1813  }
1814 
1815  // Restore current model parameter
1816  *par = current;
1817 
1818  } // endif: parameter was free
1819 
1820  } // endfor: looped over all parameters
1821 
1822  // Take square root if required
1823  if (eflux_error > 0.0) {
1824  eflux_error = std::sqrt(eflux_error);
1825  }
1826  // Return flux error
1827  return eflux_error;
1828 }
1829 
1830 
1831 /***********************************************************************//**
1832  * @brief Print model information
1833  *
1834  * @param[in] chatter Chattiness.
1835  * @return String containing model information.
1836  ***************************************************************************/
1837 std::string GModelSky::print(const GChatter& chatter) const
1838 {
1839  // Initialise result string
1840  std::string result;
1841 
1842  // Continue only if chatter is not silent
1843  if (chatter != SILENT) {
1844 
1845  // Append header
1846  result.append("=== GModelSky ===");
1847 
1848  // Determine number of parameters per type
1849  int n_spatial = (m_spatial != NULL) ? m_spatial->size() : 0;
1850  int n_spectral = (m_spectral != NULL) ? m_spectral->size() : 0;
1851  int n_temporal = (m_temporal != NULL) ? m_temporal->size() : 0;
1852 
1853  // Append attributes
1854  result.append("\n"+print_attributes());
1855 
1856  // Append model type
1857  result.append("\n"+gammalib::parformat("Model type")+type());
1858 
1859  // Append model components
1860  result.append("\n"+gammalib::parformat("Model components"));
1861  if (n_spatial > 0) {
1862  result.append("\""+spatial()->type()+"\"");
1863  if (n_spectral > 0 || n_temporal > 0) {
1864  result.append(" * ");
1865  }
1866  }
1867  if (n_spectral > 0) {
1868  result.append("\""+spectral()->type()+"\"");
1869  if (n_temporal > 0) {
1870  result.append(" * ");
1871  }
1872  }
1873  if (n_temporal > 0) {
1874  result.append("\""+temporal()->type()+"\"");
1875  }
1876 
1877  // Append parameters
1878  result.append("\n"+gammalib::parformat("Number of parameters"));
1879  result.append(gammalib::str(size()));
1880  result.append("\n"+gammalib::parformat("Number of spatial par's"));
1881  result.append(gammalib::str(n_spatial));
1882  for (int i = 0; i < n_spatial; ++i) {
1883  result.append("\n"+(*spatial())[i].print());
1884  }
1885  result.append("\n"+gammalib::parformat("Number of spectral par's"));
1886  result.append(gammalib::str(n_spectral));
1887  for (int i = 0; i < n_spectral; ++i) {
1888  result.append("\n"+(*spectral())[i].print());
1889  }
1890  result.append("\n"+gammalib::parformat("Number of temporal par's"));
1891  result.append(gammalib::str(n_temporal));
1892  for (int i = 0; i < n_temporal; ++i) {
1893  result.append("\n"+(*temporal())[i].print());
1894  }
1895  result.append("\n"+gammalib::parformat("Number of scale par's"));
1896  result.append(gammalib::str(scales()));
1897  for (int i = 0; i < scales(); ++i) {
1898  result.append("\n"+scale(i).print());
1899  }
1900 
1901  } // endif: chatter was not silent
1902 
1903  // Return result
1904  return result;
1905 }
1906 
1907 
1908 /*==========================================================================
1909  = =
1910  = Private methods =
1911  = =
1912  ==========================================================================*/
1913 
1914 /***********************************************************************//**
1915  * @brief Initialise class members
1916  ***************************************************************************/
1918 {
1919  // Initialise members
1920  m_type.clear();
1921  m_spatial = NULL;
1922  m_spectral = NULL;
1923  m_temporal = NULL;
1924 
1925  // Return
1926  return;
1927 }
1928 
1929 
1930 /***********************************************************************//**
1931  * @brief Copy class members
1932  *
1933  * @param[in] model Sky model.
1934  ***************************************************************************/
1936 {
1937  // Copy attributes
1938  m_type = model.m_type;
1939 
1940  // Clone model components
1941  m_spatial = (model.m_spatial != NULL) ? model.m_spatial->clone() : NULL;
1942  m_spectral = (model.m_spectral != NULL) ? model.m_spectral->clone() : NULL;
1943  m_temporal = (model.m_temporal != NULL) ? model.m_temporal->clone() : NULL;
1944 
1945  // Set parameter pointers
1946  set_pointers();
1947 
1948  // Return
1949  return;
1950 }
1951 
1952 
1953 /***********************************************************************//**
1954  * @brief Delete class members
1955  ***************************************************************************/
1957 {
1958  // Free memory
1959  if (m_spatial != NULL) delete m_spatial;
1960  if (m_spectral != NULL) delete m_spectral;
1961  if (m_temporal != NULL) delete m_temporal;
1962 
1963  // Signal free pointers
1964  m_spatial = NULL;
1965  m_spectral = NULL;
1966  m_temporal = NULL;
1967 
1968  // Return
1969  return;
1970 }
1971 
1972 
1973 /***********************************************************************//**
1974  * @brief Set parameter pointers
1975  *
1976  * Gathers all parameter pointers from the model into a flat array of model
1977  * pointers.
1978  ***************************************************************************/
1980 {
1981  // Clear parameters
1982  m_pars.clear();
1983 
1984  // Determine the number of parameters
1985  int n_spatial = (spatial() != NULL) ? spatial()->size() : 0;
1986  int n_spectral = (spectral() != NULL) ? spectral()->size() : 0;
1987  int n_temporal = (temporal() != NULL) ? temporal()->size() : 0;
1988  int n_pars = n_spatial + n_spectral + n_temporal;
1989  int n_scales = m_scales.size();
1990 
1991  // Continue only if there are parameters
1992  if (n_pars > 0) {
1993 
1994  // Gather spatial parameter pointers
1995  for (int i = 0; i < n_spatial; ++i) {
1996  m_pars.push_back(&((*spatial())[i]));
1997  }
1998 
1999  // Gather spectral parameters
2000  for (int i = 0; i < n_spectral; ++i) {
2001  m_pars.push_back(&((*spectral())[i]));
2002  }
2003 
2004  // Gather temporal parameters
2005  for (int i = 0; i < n_temporal; ++i) {
2006  m_pars.push_back(&((*temporal())[i]));
2007  }
2008 
2009  }
2010 
2011  // Append scales if they exist
2012  if (n_scales > 0) {
2013  for (int i = 0; i < n_scales; ++i) {
2014  m_pars.push_back(&m_scales[i]);
2015  }
2016  }
2017 
2018  // Return
2019  return;
2020 }
2021 
2022 
2023 /***********************************************************************//**
2024  * @brief Set model type based on spatial model component
2025  *
2026  * Set the model type depending on the class type of the spatial model
2027  * component.
2028  *
2029  * @todo A method could be implemented in the GModelSpatial class that
2030  * determines the model type. This is however not very critical.
2031  ***************************************************************************/
2033 {
2034  // Clear model type
2035  m_type.clear();
2036 
2037  // Continue only if we have a spatial model component
2038  if (m_spatial != NULL) {
2039 
2040  // Is spatial model a point source?
2041  if (dynamic_cast<const GModelSpatialPointSource*>(m_spatial) != NULL) {
2042  m_type = "PointSource";
2043  }
2044 
2045  // ... otherwise, is spatial model a radial source?
2046  else if (dynamic_cast<const GModelSpatialRadial*>(m_spatial) != NULL) {
2047  m_type = "ExtendedSource";
2048  }
2049 
2050  // ... otherwise, is spatial model an elliptical source?
2051  else if (dynamic_cast<const GModelSpatialElliptical*>(m_spatial) != NULL) {
2052  m_type = "ExtendedSource";
2053  }
2054 
2055  // ... otherwise, is spatial model a composite source?
2056  else if (dynamic_cast<const GModelSpatialComposite*>(m_spatial) != NULL) {
2057  m_type = "CompositeSource";
2058  }
2059 
2060  // ... otherwise we have a diffuse model
2061  else {
2062  m_type = "DiffuseSource";
2063  }
2064 
2065  } // endif: there was a spatial model component
2066 
2067  // Return
2068  return;
2069 }
2070 
2071 
2072 /***********************************************************************//**
2073  * @brief Signal all parameters that have analytical gradients for a given
2074  * observation
2075  *
2076  * @param[in] obs Observation.
2077  *
2078  * This method signals all spectral, temporal and model scaling parameters
2079  * for which analytical gradients exist for a given observation. Analytical
2080  * gradients will be signalled for all of these parameters that are free
2081  * and for which `has_grad()` is `true`.
2082  *
2083  * Note that spatial analytical gradients are excluded. These will be
2084  * signalled by the underlying convolution method as analytical gradients
2085  * are not (yet) supported by all convolution methods.
2086  ***************************************************************************/
2088 {
2089  // If the model has a spectral component then signal all free parameters
2090  // with analytical gradients
2091  if (spectral() != NULL) {
2092  for (int i = 0; i < spectral()->size(); ++i) {
2093  GModelPar& par = (*(spectral()))[i];
2094  if (par.is_free() && par.has_grad()) {
2095  obs.computed_gradient(*this, par);
2096  }
2097  }
2098  }
2099 
2100  // If the model has a temporal component then signal all free parameters
2101  // with analytical gradients
2102  if (temporal() != NULL) {
2103  for (int i = 0; i < temporal()->size(); ++i) {
2104  GModelPar& par = (*(temporal()))[i];
2105  if (par.is_free() && par.has_grad()) {
2106  obs.computed_gradient(*this, par);
2107  }
2108  }
2109  }
2110 
2111  // If the model has a scales then signal all free parameters with
2112  // analytical gradients
2113  if (has_scales()) {
2114  for (int i = 0; i < scales(); ++i) {
2115  const GModelPar& par = scale(i);
2116  if (par.name() == obs.instrument()) {
2117  if (par.is_free() && par.has_grad()) {
2118  obs.computed_gradient(*this, par);
2119  }
2120  }
2121  }
2122  }
2123 
2124  // Return
2125  return;
2126 }
2127 
2128 
2129 /***********************************************************************//**
2130  * @brief Return pointer to spatial model from XML element
2131  *
2132  * @param[in] spatial XML element.
2133  * @return Pointer to spatial model.
2134  *
2135  * Returns pointer to spatial model that is defined in an XML element.
2136  ***************************************************************************/
2138 {
2139  // Get spatial model
2140  GModelSpatialRegistry registry;
2141  GModelSpatial* ptr = registry.alloc(spatial);
2142 
2143  // Return pointer
2144  return ptr;
2145 }
2146 
2147 
2148 /***********************************************************************//**
2149  * @brief Return pointer to spectral model from XML element
2150  *
2151  * @param[in] spectral XML element.
2152  * @return Pointer to spectral model.
2153  *
2154  * Returns pointer to spectral model that is defined in an XML element.
2155  ***************************************************************************/
2157 {
2158  // Get spectral model
2159  GModelSpectralRegistry registry;
2160  GModelSpectral* ptr = registry.alloc(spectral);
2161 
2162  // Return pointer
2163  return ptr;
2164 }
2165 
2166 
2167 /***********************************************************************//**
2168  * @brief Return pointer to temporal model from XML element
2169  *
2170  * @param[in] temporal XML element.
2171  * @return Pointer to temporal model.
2172  *
2173  * Returns pointer to temporal model that is defined in an XML element.
2174  ***************************************************************************/
2176 {
2177  // Get temporal model
2178  GModelTemporalRegistry registry;
2179  GModelTemporal* ptr = registry.alloc(temporal);
2180 
2181  // Return pointer
2182  return ptr;
2183 }
2184 
2185 
2186 /***********************************************************************//**
2187  * @brief Verifies if model has all components
2188  *
2189  * @return True is model is valid, false otherwise.
2190  ***************************************************************************/
2191 bool GModelSky::valid_model(void) const
2192 {
2193  // Set result
2194  bool result = ((m_spatial != NULL) &&
2195  (m_spectral != NULL) &&
2196  (m_temporal != NULL));
2197 
2198  // Return result
2199  return result;
2200 }
2201 
2202 
2203 /***********************************************************************//**
2204  * @brief Integration kernel for flux_kern() class
2205  *
2206  * @param[in] x Function value.
2207  * @return Photon flux (erg/cm2/s/MeV)
2208  *
2209  * This method implements the integration kernel needed for the flux()
2210  * method.
2211  ***************************************************************************/
2212 double GModelSky::flux_kern::eval(const double& x)
2213 {
2214  // Set energy
2215  GEnergy eng;
2216  double expx = std::exp(x);
2217  eng.MeV(expx);
2218 
2219  // Get photon flux within sky region -> (ph/cm2/s)
2220  double value = m_parent->spatial()->flux(*m_region, eng);
2221 
2222  // Multiply by spectral model -> (ph/cm2/s/MeV)
2223  value *= m_parent->spectral()->eval(eng);
2224 
2225  // Correct for variable substitution
2226  value *= expx;
2227 
2228  // Return value
2229  return value;
2230 }
2231 
2232 
2233 /***********************************************************************//**
2234  * @brief Integration kernel for eflux_kern() class
2235  *
2236  * @param[in] x Function value.
2237  * @return Energy flux (erg/cm2/s/MeV)
2238  *
2239  * This method implements the integration kernel needed for the eflux()
2240  * method.
2241  ***************************************************************************/
2242 double GModelSky::eflux_kern::eval(const double& x)
2243 {
2244  // Set energy
2245  GEnergy eng;
2246  double expx = std::exp(x);
2247  eng.MeV(expx);
2248 
2249  // Get photon flux within sky region -> (ph/cm2/s)
2250  double value = m_parent->spatial()->flux(*m_region, eng);
2251 
2252  // Multiply by spectral model -> (ph/cm2/s/MeV)
2253  value *= m_parent->spectral()->eval(eng);
2254 
2255  // Correct for variable substitution
2256  value *= expx;
2257 
2258  // Multiply by energy in MeV -> (MeV/cm2/s/MeV)
2259  value *= expx * gammalib::MeV2erg;
2260 
2261  // Return value
2262  return value;
2263 }
virtual void write(GXmlElement &xml) const =0
Constant temporal model class interface definition.
const double & factor_error(void) const
Return parameter factor error.
Abstract model class.
Definition: GModel.hpp:100
double dist_deg(const GSkyDir &dir) const
Compute angular distance between sky directions in degrees.
Definition: GSkyDir.hpp:286
void free_members(void)
Delete class members.
Definition: GModel.cpp:672
int scales(void) const
Return number of scale parameters in model.
Definition: GModel.hpp:247
void set_type(void)
Set model type based on spatial model component.
Definition: GModelSky.cpp:2032
void computed_gradient(const GModel &model, const GModelPar &par) const
Signals that an analytical gradient was computed for a model parameter.
virtual GSkyDir mc(const GEnergy &energy, const GTime &time, GRan &ran) const =0
GModelTemporal * m_temporal
Temporal model.
Definition: GModelSky.hpp:239
virtual GTimes mc(const double &rate, const GTime &tmin, const GTime &tmax, GRan &ran) const =0
virtual GEnergy mc(const GEnergy &emin, const GEnergy &emax, const GTime &time, GRan &ran) const =0
double norm(const GVector &vector)
Computes vector norm.
Definition: GVector.cpp:864
double romberg(std::vector< double > bounds, const int &order=5)
Perform Romberg integration.
Definition: GIntegral.cpp:381
Point source spatial model class interface definition.
int size(void) const
Return number of times.
Definition: GTimes.hpp:100
virtual double eval(const GPhoton &photon, const bool &gradients=false) const =0
const std::string & name(void) const
Return parameter name.
std::string print_attributes(void) const
Print model attributes.
Definition: GModel.cpp:770
Abstract spectral model base class.
Sparse matrix class interface definition.
void copy_members(const GModelSky &model)
Copy class members.
Definition: GModelSky.cpp:1935
Interface definition for the model registry class.
virtual std::string instrument(void) const =0
GModelSpectral * spectral(void) const
Return spectral model component.
Definition: GModelSky.hpp:311
GModelTemporal * xml_temporal(const GXmlElement &temporal) const
Return pointer to temporal model from XML element.
Definition: GModelSky.cpp:2175
double eflux_error(const GEnergy &emin, const GEnergy &emax) const
Return sky model energy flux error.
Definition: GModelSky.cpp:1575
int size(void) const
Return number of parameters.
Abstract temporal model base class.
GPhotons mc(const double &area, const GSkyDir &dir, const double &radius, const GEnergy &emin, const GEnergy &emax, const GTime &tmin, const GTime &tmax, GRan &ran) const
Return simulated list of photons.
Definition: GModelSky.cpp:882
int size(void) const
Return number of parameters.
virtual void clear(void)
Clear sky model.
Definition: GModelSky.cpp:379
virtual double eval(const GEnergy &srcEng, const GTime &srcTime=GTime(), const bool &gradients=false) const =0
GModelSpatial * xml_spatial(const GXmlElement &spatial) const
Return pointer to spatial model from XML element.
Definition: GModelSky.cpp:2137
GModelTemporal * temporal(void) const
Return temporal model component.
Definition: GModelSky.hpp:327
Abstract radial spatial model base class interface definition.
void signal_analytical_gradients(const GObservation &obs) const
Signal all parameters that have analytical gradients for a given observation.
Definition: GModelSky.cpp:2087
Abstract interface for the event classes.
Definition: GEvent.hpp:71
void set_pointers(void)
Set parameter pointers.
Definition: GModelSky.cpp:1979
XML element node class.
Definition: GXmlElement.hpp:48
Spectral model registry class definition.
virtual ~GModelSky(void)
Destructor.
Definition: GModelSky.cpp:322
virtual void write(GXmlElement &xml) const =0
Random number generator class.
Definition: GRan.hpp:44
virtual double flux(const GEnergy &emin, const GEnergy &emax) const
Returns model photon flux between emin, emax
virtual double eflux(const GEnergy &emin, const GEnergy &emax) const
Returns model energy flux between emin, emax
double MeV(void) const
Return energy in MeV.
Definition: GEnergy.cpp:321
Photon container class.
Definition: GPhotons.hpp:45
bool valid_model(void) const
Verifies if model has all components.
Definition: GModelSky.cpp:2191
virtual int elements(void) const
Return number of GXMLElement children of node.
Definition: GXmlNode.cpp:586
Time class.
Definition: GTime.hpp:55
Gammalib tools definition.
GIntegral class interface definition.
Definition: GIntegral.hpp:46
virtual GModelTemporalConst * clone(void) const
Clone constant temporal model.
Interface for the circular sky region class.
Interface definition for the spatial model registry class.
Source class definition.
const GModelSky * m_parent
Sky model.
Definition: GModelSky.hpp:220
bool is_free(void) const
Signal if parameter is free.
std::string m_type
Model type.
Definition: GModelSky.hpp:236
GModelSpectral * m_spectral
Spectral model.
Definition: GModelSky.hpp:238
std::vector< GModelPar > m_scales
Model instrument scale factors.
Definition: GModel.hpp:176
Spectral nodes model class.
double flux_error(const GEnergy &emin, const GEnergy &emax) const
Return sky model photon flux error.
Definition: GModelSky.cpp:1299
Time container class.
Definition: GTimes.hpp:45
bool is_notanumber(const double &x)
Signal if argument is not a number.
Definition: GTools.hpp:201
std::vector< GModelPar * > m_pars
Pointers to all model parameters.
Definition: GModel.hpp:178
virtual void write(GXmlElement &xml) const
Write model into XML element.
Definition: GModelSky.cpp:793
int size(void) const
Return number of parameters in model.
Definition: GModel.hpp:233
Class that handles photons.
Definition: GPhoton.hpp:47
bool is_infinite(const double &x)
Signal if argument is infinite.
Definition: GTools.hpp:184
void init_members(void)
Initialise class members.
Definition: GModel.cpp:622
Model parameter class.
Definition: GModelPar.hpp:87
void mc_cone(const GSkyRegionCircle &cone) const
Set Monte Carlo simulation cone.
GVector sqrt(const GVector &vector)
Computes square root of vector elements.
Definition: GVector.cpp:1358
const GModelSpectralNodes & spectrum(void) const
Get map cube spectrum.
bool has_grad(void) const
Signal if parameter gradient is computed analytically.
Abstract interface for the sky region class.
Definition: GSkyRegion.hpp:57
double eval(const double &x)
Integration kernel for flux_kern() class.
Definition: GModelSky.cpp:2212
const GXmlAttribute * attribute(const int &index) const
Return attribute.
const std::string & name(void) const
Return parameter name.
Definition: GModel.hpp:261
const double & factor_max(void) const
Return parameter maximum factor boundary.
GVector gradients(const GPhoton &photon)
Return parameter gradients of sky model for a given photon.
Definition: GModelSky.cpp:529
double eflux(const GEnergy &emin, const GEnergy &emax) const
Return sky model energy flux.
Definition: GModelSky.cpp:1173
Energy container class definition.
double eval(const double &x)
Integration kernel for eflux_kern() class.
Definition: GModelSky.cpp:2242
int size(void) const
Return number of parameters.
virtual double npred(const GEnergy &obsEng, const GTime &obsTime, const GObservation &obs) const
Return spatially integrated sky model.
Definition: GModelSky.cpp:648
const GSkyRegion * m_region
Sky region.
Definition: GModelSky.hpp:221
GModelTemporal * alloc(const GXmlElement &xml) const
Allocate temporal model that is found in XML element.
void autoscale(void)
Autoscale parameter.
GVector log(const GVector &vector)
Computes natural logarithm of vector elements.
Definition: GVector.cpp:1274
Spatial model registry class definition.
virtual void read(const GXmlElement &xml)
Read sky model from XML element.
Definition: GModelSky.cpp:724
GModelSpectral * alloc(const GXmlElement &xml) const
Allocate spectral model that is found in XML element.
const GTime & time(void) const
Return photon time.
Definition: GPhoton.hpp:134
GChatter
Definition: GTypemaps.hpp:33
GEnergy energy(const int &index) const
Return node energy.
GModelSpatial * alloc(const GXmlElement &xml) const
Allocate spatial model that is found in XML element.
virtual void response(const GResponse &rsp)=0
virtual double eval(const GTime &srcTime, const bool &gradients=false) const =0
bool has_max(void) const
Signal if parameter has maximum boundary.
Spatial composite model class interface definition.
Abstract observation base class.
virtual GModelSky * clone(void) const
Clone sky model.
Definition: GModelSky.cpp:399
bool has_scales(void) const
Signals that model has scales.
Definition: GModel.hpp:355
const GModelSky g_extendedsource_seed("ExtendedSource")
Definition: GModelSky.cpp:55
GModelPar & scale(const int &index)
Returns reference to scale parameter by index.
Definition: GModel.cpp:369
Interface definition for the spectral model registry class.
Interface definition for the temporal model registry class.
void write_attributes(GXmlElement &xml) const
Write model attributes.
Definition: GModel.cpp:723
virtual GModelSpectral * clone(void) const =0
Clones object.
virtual GModelSky & operator=(const GModelSky &model)
Assignment operator.
Definition: GModelSky.cpp:344
Abstract elliptical spatial model base class interface definition.
Abstract response base class definition.
double intensity(const int &index) const
Return node intensity.
const double & factor_min(void) const
Return parameter minimum factor boundary.
double value(const GPhoton &photon)
Return value of sky model for a given photon.
Definition: GModelSky.cpp:503
const GModelSky g_diffusesource_seed("DiffuseSource")
Definition: GModelSky.cpp:56
virtual GXmlElement * element(const int &index)
Return pointer to GXMLElement child.
Definition: GXmlNode.cpp:640
void reserve(const int &num)
Reserve memory for photons in container.
Definition: GPhotons.cpp:334
Temporal model registry class definition.
Sky model class interface definition.
Sky model class.
Definition: GModelSky.hpp:122
const double MeV2erg
Definition: GTools.hpp:45
Spatial map cube model class interface definition.
int nodes(void) const
Return number of nodes.
bool has_min(void) const
Signal if parameter has minimum boundary.
void init_members(void)
Initialise class members.
Definition: GModelSky.cpp:1917
virtual std::string classname(void) const =0
Return class name.
virtual double mc_norm(const GSkyDir &dir, const double &radius) const =0
void append(const GPhoton &photon)
Append photon to container.
Definition: GPhotons.cpp:219
virtual std::string type(void) const
Return sky model type.
Definition: GModelSky.hpp:264
virtual double flux(const GSkyRegion &region, const GEnergy &srcEng=GEnergy(), const GTime &srcTime=GTime()) const
Returns model flux integrated in circular sky region.
GModelSpatial * m_spatial
Spatial model.
Definition: GModelSky.hpp:237
const GEnergy & energy(void) const
Return photon energy.
Definition: GPhoton.hpp:122
Exception handler interface definition.
GModelSpatial * spatial(void) const
Return spatial model component.
Definition: GModelSky.hpp:295
GVector exp(const GVector &vector)
Computes exponential of vector elements.
Definition: GVector.cpp:1232
virtual GModelTemporal * clone(void) const =0
Clones object.
void read_attributes(const GXmlElement &xml)
Read model attributes.
Definition: GModel.cpp:684
Abstract spatial model base class.
virtual void write(GXmlElement &xml) const =0
GModelSky(void)
Void constructor.
Definition: GModelSky.cpp:91
double norm(void) const
Return normalization factor.
virtual GModel & operator=(const GModel &model)
Assignment operator.
Definition: GModel.cpp:132
Vector class.
Definition: GVector.hpp:46
virtual double flux(const GEnergy &emin, const GEnergy &emax) const =0
Model registry class definition.
virtual double eflux(const GEnergy &emin, const GEnergy &emax) const =0
virtual GXmlNode * append(const GXmlNode &node)
Append XML child node.
Definition: GXmlNode.cpp:287
std::string parformat(const std::string &s, const int &indent=0)
Convert string in parameter format.
Definition: GTools.cpp:1143
double flux(const GEnergy &emin, const GEnergy &emax) const
Return sky model photon flux.
Definition: GModelSky.cpp:1061
virtual std::string print(const GChatter &chatter=NORMAL) const
Print model information.
Definition: GModelSky.cpp:1837
Integration class interface definition.
Sky direction class.
Definition: GSkyDir.hpp:62
const GModelSky g_compositesource_seed("CompositeSource")
Definition: GModelSky.cpp:57
virtual GModelSpatial * clone(void) const =0
Clones object.
const GSkyDir & dir(void) const
Return photon sky direction.
Definition: GPhoton.hpp:110
virtual double eval(const GEvent &event, const GObservation &obs, const bool &gradients=false) const
Evaluate sky model for a given event of an observation.
Definition: GModelSky.cpp:565
const GModelSky g_pointsource_seed("PointSource")
Definition: GModelSky.cpp:54
Constant temporal model class.
void free_members(void)
Delete class members.
Definition: GModelSky.cpp:1956
const double & factor_value(void) const
Return parameter factor value.
Class that handles energies in a unit independent way.
Definition: GEnergy.hpp:48
GModelSpectral * xml_spectral(const GXmlElement &spectral) const
Return pointer to spectral model from XML element.
Definition: GModelSky.cpp:2156
std::string str(const unsigned short int &value)
Convert unsigned short integer value into string.
Definition: GTools.cpp:489