GammaLib  1.7.0.dev
 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-2019 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 "GModelRegistry.hpp"
35 #include "GModelSky.hpp"
38 #include "GModelSpatialRadial.hpp"
44 #include "GModelTemporalConst.hpp"
45 #include "GSource.hpp"
46 #include "GResponse.hpp"
47 
48 /* __ Globals ____________________________________________________________ */
49 const GModelSky g_pointsource_seed("PointSource");
50 const GModelSky g_extendedsource_seed("ExtendedSource");
51 const GModelSky g_diffusesource_seed("DiffuseSource");
52 const GModelSky g_compositesource_seed("CompositeSource");
53 const GModelRegistry g_pointsource_registry(&g_pointsource_seed);
54 const GModelRegistry g_extendedsource_registry(&g_extendedsource_seed);
55 const GModelRegistry g_diffusesource_registry(&g_diffusesource_seed);
56 const GModelRegistry g_compositesource_registry(&g_compositesource_seed);
57 
58 /* __ Method name definitions ____________________________________________ */
59 #define G_NPRED "GModelSky::npred(GEnergy&, GTime&, GObservation&)"
60 #define G_XML_SPATIAL "GModelSky::xml_spatial(GXmlElement&)"
61 #define G_XML_SPECTRAL "GModelSky::xml_spectral(GXmlElement&)"
62 #define G_XML_TEMPORAL "GModelSky::xml_temporal(GXmlElement&)"
63 #define G_INTEGRATE_TIME "GModelSky::integrate_time(GEvent&, GObservation&,"\
64  " bool)"
65 #define G_INTEGRATE_ENERGY "GModelSky::integrate_energy(GEvent&, GTime&," \
66  " GObservation&, bool)"
67 #define G_INTEGRATE_DIR "GModelSky::integrate_dir(GEvent&, GEnergy&, GTime&,"\
68  " GObservation, bool)"
69 
70 /* __ Macros _____________________________________________________________ */
71 
72 /* __ Coding definitions _________________________________________________ */
73 
74 /* __ Debug definitions __________________________________________________ */
75 //#define G_DUMP_MC //!< Dump MC information
76 //#define G_DUMP_MC_DETAIL //!< Dump detailed MC information
77 
78 
79 /*==========================================================================
80  = =
81  = Constructors/destructors =
82  = =
83  ==========================================================================*/
84 
85 /***********************************************************************//**
86  * @brief Void constructor
87  *
88  * Construct an empty sky model. An empty sky model has no type.
89  ***************************************************************************/
91 {
92  // Initialise members
93  init_members();
94 
95  // Return
96  return;
97 }
98 
99 
100 /***********************************************************************//**
101  * @brief Type constructor
102  *
103  * @param[in] type Model type.
104  *
105  * Construct an empty sky model of the specified @p type. This constructor
106  * does basically the same than the void constructor with the only difference
107  * that the m_type member of the class is set to the specified @p type
108  * string.
109  ***************************************************************************/
110 GModelSky::GModelSky(const std::string& type) : GModel()
111 {
112  // Initialise members
113  init_members();
114 
115  // Set model type
116  m_type = type;
117 
118  // Return
119  return;
120 }
121 
122 
123 /***********************************************************************//**
124  * @brief XML constructor
125  *
126  * @param[in] xml XML element.
127  *
128  * Constructs a sky model from the model information that is found in the
129  * @p xml element. Only the spatial and spectral information is handled by
130  * this constructor. For the temporal component, a constant model of type
131  * GModelTemporalConst will be allocated.
132  ***************************************************************************/
134 {
135  // Initialise members
136  init_members();
137 
138  // Get pointers on spectrum and spatial model
139  const GXmlElement* spec = xml.element("spectrum", 0);
140  const GXmlElement* spat = xml.element("spatialModel", 0);
141 
142  // Allocate constant
144 
145  // Clone spatial and spectral models
146  m_spatial = xml_spatial(*spat);
147  m_spectral = xml_spectral(*spec);
148  m_temporal = temporal.clone();
149 
150  // Set parameter pointers
151  set_pointers();
152 
153  // Set model type dependent on spatial model type
154  set_type();
155 
156  // Return
157  return;
158 }
159 
160 
161 /***********************************************************************//**
162  * @brief Construct sky model from spatial and spectral XML elements
163  *
164  * @param[in] spatial Spatial XML element.
165  * @param[in] spectral Spectral XML element.
166  *
167  * Constructs a sky model from the @p spatial and @p spectral information
168  * that is found in the respective XML elements. For the temporal component,
169  * a constant model of type GModelTemporalConst will be allocated.
170  ***************************************************************************/
171 GModelSky::GModelSky(const GXmlElement& spatial, const GXmlElement& spectral) :
172  GModel()
173 {
174  // Initialise private members for clean destruction
175  init_members();
176 
177  // Allocate constant
179 
180  // Clone spatial and spectral models
181  m_spatial = xml_spatial(spatial);
182  m_spectral = xml_spectral(spectral);
183  m_temporal = temporal.clone();
184 
185  // Set parameter pointers
186  set_pointers();
187 
188  // Set model type dependent on spatial model type
189  set_type();
190 
191  // Return
192  return;
193 }
194 
195 
196 /***********************************************************************//**
197  * @brief Construct sky model from spatial, spectral and temporal XML elements
198  *
199  * @param[in] spatial Spatial XML element.
200  * @param[in] spectral Spectral XML element.
201  * @param[in] temporal Temporal XML element.
202  *
203  * Constructs a sky model from the @p spatial, @p spectral and @p temporal
204  * information that is found in the respective XML elements.
205  ***************************************************************************/
207  const GXmlElement& spectral,
208  const GXmlElement& temporal) :
209  GModel()
210 {
211  // Initialise private members for clean destruction
212  init_members();
213 
214  // Clone spatial and spectral models
215  m_spatial = xml_spatial(spatial);
216  m_spectral = xml_spectral(spectral);
217  m_temporal = xml_temporal(temporal);
218 
219  // Set parameter pointers
220  set_pointers();
221 
222  // Set model type dependent on spatial model type
223  set_type();
224 
225  // Return
226  return;
227 }
228 
229 
230 /***********************************************************************//**
231  * @brief Construct sky model from spatial and spectral components
232  *
233  * @param[in] spatial Spatial model component.
234  * @param[in] spectral Spectral model component.
235  *
236  * Constructs a sky model from @p spatial and @p spectral model components.
237  * For the temporal component, a constant model of type GModelTemporalConst
238  * will be allocated.
239  ***************************************************************************/
241  const GModelSpectral& spectral) :
242  GModel()
243 {
244  // Initialise members
245  init_members();
246 
247  // Allocate temporal constant model
249 
250  // Clone model components
251  m_spatial = spatial.clone();
252  m_spectral = spectral.clone();
253  m_temporal = temporal.clone();
254 
255  // Set parameter pointers
256  set_pointers();
257 
258  // Set model type dependent on spatial model type
259  set_type();
260 
261  // Return
262  return;
263 }
264 
265 
266 /***********************************************************************//**
267  * @brief Construct sky model from spatial, spectral and temporal components
268  *
269  * @param[in] spatial Spatial model component.
270  * @param[in] spectral Spectral model component.
271  * @param[in] temporal Temporal model component.
272  *
273  * Constructs a sky model from @p spatial, @p spectral and @p temporal model
274  * components.
275  ***************************************************************************/
277  const GModelSpectral& spectral,
278  const GModelTemporal& temporal) :
279  GModel()
280 {
281  // Initialise members
282  init_members();
283 
284  // Clone model components
285  m_spatial = spatial.clone();
286  m_spectral = spectral.clone();
287  m_temporal = temporal.clone();
288 
289  // Set parameter pointers
290  set_pointers();
291 
292  // Set model type dependent on spatial model type
293  set_type();
294 
295  // Return
296  return;
297 }
298 
299 
300 /***********************************************************************//**
301  * @brief Copy constructor
302  *
303  * @param[in] model Sky model.
304  ***************************************************************************/
305 GModelSky::GModelSky(const GModelSky& model) : GModel(model)
306 {
307  // Initialise private members for clean destruction
308  init_members();
309 
310  // Copy members
311  copy_members(model);
312 
313  // Return
314  return;
315 }
316 
317 
318 /***********************************************************************//**
319  * @brief Destructor
320  ***************************************************************************/
322 {
323  // Free members
324  free_members();
325 
326  // Return
327  return;
328 }
329 
330 
331 /*==========================================================================
332  = =
333  = Operators =
334  = =
335  ==========================================================================*/
336 
337 /***********************************************************************//**
338  * @brief Assignment operator
339  *
340  * @param[in] model Sky model.
341  * @return Sky model.
342  ***************************************************************************/
344 {
345  // Execute only if object is not identical
346  if (this != &model) {
347 
348  // Copy base class members
349  this->GModel::operator=(model);
350 
351  // Free members
352  free_members();
353 
354  // Initialise private members for clean destruction
355  init_members();
356 
357  // Copy members
358  copy_members(model);
359 
360  } // endif: object was not identical
361 
362  // Return
363  return *this;
364 }
365 
366 
367 /*==========================================================================
368  = =
369  = Public methods =
370  = =
371  ==========================================================================*/
372 
373 /***********************************************************************//**
374  * @brief Clear sky model
375  *
376  * This method properly resets the sky model to an initial state.
377  ***************************************************************************/
379 {
380  // Free class members
381  free_members();
382  this->GModel::free_members();
383 
384  // Initialise members
385  this->GModel::init_members();
386  init_members();
387 
388  // Return
389  return;
390 }
391 
392 
393 /***********************************************************************//**
394  * @brief Clone sky model
395  *
396  * @return Pointer to deep copy of sky model.
397  ***************************************************************************/
399 {
400  // Clone sky model
401  return new GModelSky(*this);
402 }
403 
404 
405 /***********************************************************************//**
406  * @brief Set spatial model component
407  *
408  * @param[in] spatial Pointer to spatial model component.
409  *
410  * Sets the spatial model component of the model.
411  ***************************************************************************/
412 void GModelSky::spatial(const GModelSpatial* spatial)
413 {
414  // Free spatial model component only if it differs from current
415  // component. This prevents unintential deallocation of the argument
416  if ((m_spatial != NULL) && (m_spatial != spatial)) {
417  delete m_spatial;
418  }
419 
420  // Clone spatial model component if it exists, otherwise set pointer
421  // to NULL
422  m_spatial = (spatial != NULL) ? spatial->clone() : NULL;
423 
424  // Set parameter pointers
425  set_pointers();
426 
427  // Set model type dependent on spatial model type
428  set_type();
429 
430  // Return
431  return;
432 }
433 
434 
435 /***********************************************************************//**
436  * @brief Set spectral model component
437  *
438  * @param[in] spectral Pointer to spectral model component.
439  *
440  * Sets the spectral model component of the model.
441  ***************************************************************************/
442 void GModelSky::spectral(const GModelSpectral* spectral)
443 {
444  // Free spectral model component only if it differs from current
445  // component. This prevents unintential deallocation of the argument
446  if ((m_spectral != NULL) && (m_spectral != spectral)) {
447  delete m_spectral;
448  }
449 
450  // Clone spectral model component if it exists, otherwise set pointer
451  // to NULL
452  m_spectral = (spectral != NULL) ? spectral->clone() : NULL;
453 
454  // Set parameter pointers
455  set_pointers();
456 
457  // Return
458  return;
459 }
460 
461 
462 /***********************************************************************//**
463  * @brief Set temporal model component
464  *
465  * @param[in] temporal Pointer to temporal model component.
466  *
467  * Sets the temporal model component of the model.
468  ***************************************************************************/
469 void GModelSky::temporal(const GModelTemporal* temporal)
470 {
471  // Free temporal model component only if it differs from current
472  // component. This prevents unintential deallocation of the argument
473  if ((m_temporal != NULL) && (m_temporal != temporal)) {
474  delete m_temporal;
475  }
476 
477  // Clone temporal model component if it exists, otherwise set pointer
478  // to NULL
479  m_temporal = (temporal != NULL) ? temporal->clone() : NULL;
480 
481  // Set parameter pointers
482  set_pointers();
483 
484  // Return
485  return;
486 }
487 
488 
489 /***********************************************************************//**
490  * @brief Return value of sky model for a given photon
491  *
492  * @param[in] photon Photon.
493  * @return Value of sky model.
494  *
495  * Returns the value of the sky model for a given @p photon. If no model
496  * components exist the model will return a value of 1.
497  *
498  * @todo We probably should return a value of 0 is no model components exist.
499  * But we may also make sure that the model never has NULL pointers, which
500  * would avoid all the pointer checks.
501  ***************************************************************************/
502 double GModelSky::value(const GPhoton& photon)
503 {
504  // Initialise source model
505  double source = 1.0;
506 
507  // Evaluate source model
508  if (m_spatial != NULL) source *= m_spatial->eval(photon);
509  if (m_spectral != NULL) source *= m_spectral->eval(photon.energy(),
510  photon.time());
511  if (m_temporal != NULL) source *= m_temporal->eval(photon.time());
512 
513  // Return
514  return source;
515 }
516 
517 
518 /***********************************************************************//**
519  * @brief Return parameter gradients of sky model for a given photon
520  *
521  * @param[in] photon Photon.
522  * @return Vector of parameter gradients
523  *
524  * Returns a vector of parameter gradients for the sky model for a given
525  * @p photon. If there are no parameters in the sky model, an empty vector
526  * will be returned.
527  ***************************************************************************/
529 {
530  // Evaluate source model gradients
531  if (m_spatial != NULL) m_spatial->eval(photon, true);
532  if (m_spectral != NULL) m_spectral->eval(photon.energy(), photon.time(),
533  true);
534  if (m_temporal != NULL) m_temporal->eval(photon.time(), true);
535 
536  // Set vector of gradients
538  if (size() > 0) {
539  gradients = GVector(size());
540  for (int i = 0; i < size(); ++i) {
541  gradients[i] = m_pars[i]->factor_gradient();
542  }
543  }
544 
545  // Return gradients
546  return gradients;
547 }
548 
549 
550 /***********************************************************************//**
551  * @brief Evaluate sky model for a given event of an observation
552  *
553  * @param[in] event Observed event.
554  * @param[in] obs Observation.
555  * @param[in] gradients Compute gradients?
556  * @return Value of sky model.
557  *
558  * Evalutes the value of the sky model for an @p event of a specific
559  * observation @p obs.
560  *
561  * If the @p gradients flag is true the method will also compute the
562  * parameter gradients for all model parameters.
563  ***************************************************************************/
564 double GModelSky::eval(const GEvent& event,
565  const GObservation& obs,
566  const bool& gradients) const
567 {
568  // Evaluate function
569  double value = obs.response()->convolve(*this, event, obs, gradients);
570 
571  // Return
572  return value;
573 }
574 
575 
576 /***********************************************************************//**
577  * @brief Return spatially integrated sky model
578  *
579  * @param[in] obsEng Measured photon energy.
580  * @param[in] obsTime Measured photon arrival time.
581  * @param[in] obs Observation.
582  * @return Spatially integrated sky model.
583  *
584  * Computes
585  *
586  * \f[
587  * N_{\rm pred}(E',t') = \int_{\rm ROI}
588  * P(p',E',t') \, dp' \, dE'
589  * \f]
590  *
591  * of the event probability
592  *
593  * \f[
594  * P(p',E',t') = \int \int \int
595  * S(p,E,t) \times R(p',E',t'|p,E,t) \, dp \, dE \, dt
596  * \f]
597  *
598  * where
599  * \f$S(p,E,t)\f$ is the source model,
600  * \f$R(p',E',t'|p,E,t)\f$ is the instrument response function,
601  * \f$p'\f$ is the measured photon direction,
602  * \f$E'\f$ is the measured photon energy,
603  * \f$t'\f$ is the measured photon arrival time,
604  * \f$p\f$ is the true photon arrival direction,
605  * \f$E\f$ is the true photon energy, and
606  * \f$t\f$ is the true photon arrival time.
607  *
608  * The method calls the GResponse::nroi() method that does the relevant
609  * integration.
610  ***************************************************************************/
611 double GModelSky::npred(const GEnergy& obsEng,
612  const GTime& obsTime,
613  const GObservation& obs) const
614 {
615  // Initialise result
616  double npred = 0.0;
617 
618  // Continue only if model is valid)
619  if (valid_model()) {
620 
621  // Compute Nroi
622  npred = obs.response()->nroi(*this, obsEng, obsTime, obs);
623 
624  // Compile option: Check for NaN/Inf
625  #if defined(G_NAN_CHECK)
626  if (gammalib::is_notanumber(npred) || gammalib::is_infinite(npred)) {
627  std::cout << "*** ERROR: GModelSky::npred:";
628  std::cout << " NaN/Inf encountered";
629  std::cout << " (npred=" << npred;
630  std::cout << ", obsEng=" << obsEng;
631  std::cout << ", obsTime=" << obsTime;
632  std::cout << ")" << std::endl;
633  }
634  #endif
635 
636  } // endif: model was valid
637 
638  // Return npred
639  return npred;
640 }
641 
642 
643 /***********************************************************************//**
644  * @brief Read sky model from XML element
645  *
646  * @param[in] xml XML element.
647  *
648  * Reads the sky model from an XML element. The XML element is expected to
649  * respect the following format:
650  *
651  * <source name=".." type=".." instrument=".." id=".." ts="..">
652  * <spectrum type="..">
653  * ..
654  * </spectrum>
655  * <spatialModel type="..">
656  * ..
657  * </spatialModel>
658  * <temporal type="..">
659  * ..
660  * </temporal>
661  * </source>
662  *
663  * The temporal element is optional. In no temporal element is specified a
664  * constant component with unity normalization will be assumed.
665  *
666  * Optionally, the model may also contain scale parameters following the
667  * format:
668  *
669  * <source name=".." type=".." instrument=".." id=".." ts="..">
670  * <spectrum type="..">
671  * ..
672  * </spectrum>
673  * <spatialModel type="..">
674  * ..
675  * </spatialModel>
676  * <temporal type="..">
677  * ..
678  * </temporal>
679  * <scaling>
680  * <instrument name=".." scale="1.0" min="0.1" max="10.0" value="1.0" free="0"/>
681  * <instrument name=".." scale="1.0" min="0.1" max="10.0" value="0.5" free="0"/>
682  * </scaling>
683  * </source>
684  *
685  * (see GModel::read_scales() for more information on instrument scales).
686  ***************************************************************************/
687 void GModelSky::read(const GXmlElement& xml)
688 {
689  // Clear sky model
690  clear();
691 
692  // Get pointers on spectrum and spatial model
693  const GXmlElement* spec = xml.element("spectrum", 0);
694  const GXmlElement* spat = xml.element("spatialModel", 0);
695 
696  // Set spatial and spectral models
697  m_spatial = xml_spatial(*spat);
698  m_spectral = xml_spectral(*spec);
699 
700  // Handle optional temporal model
701  if (xml.elements("temporal") > 0) {
702  const GXmlElement* temp = xml.element("temporal", 0);
703  m_temporal = xml_temporal(*temp);
704  }
705  else {
707  m_temporal = temporal.clone();
708  }
709 
710  // Read model attributes
711  read_attributes(xml);
712 
713  // Set parameter pointers
714  set_pointers();
715 
716  // Set model type dependent on spatial model type
717  set_type();
718 
719  // Return
720  return;
721 }
722 
723 
724 /***********************************************************************//**
725  * @brief Write model into XML element
726  *
727  * @param[in] xml Source library.
728  *
729  * Writes the sky model into an XML source library. The format written to
730  * the @p xml element is as follows:
731  *
732  * <source name=".." type=".." instrument=".." id=".." ts="..">
733  * <spectrum type="..">
734  * ..
735  * </spectrum>
736  * <spatialModel type="..">
737  * ..
738  * </spatialModel>
739  * <temporal type="..">
740  * ..
741  * </temporal>
742  * <scaling>
743  * <instrument name=".." scale="1.0" min="0.1" max="10.0" value="1.0" free="0"/>
744  * <instrument name=".." scale="1.0" min="0.1" max="10.0" value="0.5" free="0"/>
745  * </scaling>
746  * </source>
747  *
748  * For compatibility reasons the temporal element will only be written if it
749  * is a non-constant component or a constant component with a normalization
750  * that differs from unity.
751  *
752  * The scaling element will only be written optionally in case that instrument
753  * dependent scaling factors exist (see GModel::write_scales() for more
754  * information on instrument scales).
755  ***************************************************************************/
757 {
758  // Initialise pointer on source
759  GXmlElement* src = NULL;
760 
761  // Search corresponding source
762  int n = xml.elements("source");
763  for (int k = 0; k < n; ++k) {
764  GXmlElement* element = static_cast<GXmlElement*>(xml.element("source", k));
765  if (element->attribute("name") == name()) {
766  src = element;
767  break;
768  }
769  }
770 
771  // If the temporal model is not a constant with unit normalization then
772  // set cons to a NULL pointer
773  GModelTemporalConst* cons = dynamic_cast<GModelTemporalConst*>(temporal());
774  if (cons != NULL) {
775  if (cons->norm() != 1.0) {
776  cons = NULL;
777  }
778  }
779 
780  // If no source with corresponding name was found then append one
781  if (src == NULL) {
782  src = xml.append("source");
783  if (spectral() != NULL) src->append(GXmlElement("spectrum"));
784  if (spatial() != NULL) src->append(GXmlElement("spatialModel"));
785  if (temporal() != NULL && cons == NULL) src->append(GXmlElement("temporal"));
786  }
787 
788  // Write spectral model
789  if (spectral() != NULL) {
790  GXmlElement* spec = src->element("spectrum", 0);
791  spectral()->write(*spec);
792  }
793 
794  // Write spatial model
795  if (spatial() != NULL) {
796  GXmlElement* spat = src->element("spatialModel", 0);
797  spatial()->write(*spat);
798  }
799 
800  // Write temporal model (only if not a constant with unit normalization
801  // factor)
802  if (temporal() != NULL && cons == NULL) {
803  GXmlElement* temp = src->element("temporal", 0);
804  temporal()->write(*temp);
805  }
806 
807  // Write model attributes
808  write_attributes(*src);
809 
810  // Return
811  return;
812 }
813 
814 
815 /***********************************************************************//**
816  * @brief Return simulated list of photons
817  *
818  * @param[in] area Simulation surface area (cm2).
819  * @param[in] dir Centre of simulation cone.
820  * @param[in] radius Radius of simulation cone (deg).
821  * @param[in] emin Minimum photon energy.
822  * @param[in] emax Maximum photon energy.
823  * @param[in] tmin Minimum photon arrival time.
824  * @param[in] tmax Maximum photon arrival time.
825  * @param[in,out] ran Random number generator.
826  * @return List of photons
827  *
828  * Returns a list of photons that has been derived by Monte Carlo simulation
829  * from the model. A simulation region is define by specification of
830  * - a simulation cone, which is a circular region on the sky defined by
831  * a centre direction @p dir and a @p radius,
832  * - an energy range [@p emin, @p emax], and
833  * - a time interval [@p tmin, @p tmax].
834  *
835  * Only photons with parameters in the simulation region will be returned
836  * by the method.
837  *
838  * The simulation cone may eventually cover the entire sky (by setting
839  * the radius to 180 degrees), yet simulations will be more efficient if
840  * only the sky region will be simulated that is actually observed by the
841  * telescope.
842  *
843  * @todo Implement unique model ID to assign as Monte Carlo ID
844  ***************************************************************************/
845 GPhotons GModelSky::mc(const double& area,
846  const GSkyDir& dir, const double& radius,
847  const GEnergy& emin, const GEnergy& emax,
848  const GTime& tmin, const GTime& tmax,
849  GRan& ran) const
850 {
851  // Allocate photons
852  GPhotons photons;
853 
854  // Continue only if model is valid)
855  if (valid_model()) {
856 
857  // Determine the spatial model normalization within the simulation
858  // cone and check whether the model will produce any photons in that
859  // cone.
860  double norm = m_spatial->mc_norm(dir, radius);
861  bool use_model = (norm > 0.0) ? true : false;
862 
863  // Continue only if model overlaps with simulation region
864  if (use_model) {
865 
866  // Initialise de-allocation flag
867  bool free_spectral = false;
868 
869  // Set pointer to spectral model
871 
872  // If the spectral model is a diffuse cube then create a node
873  // function spectral model that is the product of the diffuse
874  // cube node function and the spectral model evaluated at the
875  // energies of the node function
877  dynamic_cast<GModelSpatialDiffuseCube*>(m_spatial);
878  if (cube != NULL) {
879 
880  // Set MC cone
881  cube->set_mc_cone(dir, radius);
882 
883  // Allocate node function to replace the spectral component
884  GModelSpectralNodes* nodes =
885  new GModelSpectralNodes(cube->spectrum());
886  for (int i = 0; i < nodes->nodes(); ++i) {
887  GEnergy energy = nodes->energy(i);
888  double intensity = nodes->intensity(i);
889  double value = m_spectral->eval(energy, tmin);
890  nodes->intensity(i, value*intensity);
891  }
892 
893  // Signal that node function needs to be de-allocated later
894  free_spectral = true;
895 
896  // Set the spectral model pointer to the node function
897  spectral = nodes;
898 
899  // Kluge: if there are no nodes then the spectral->flux method
900  // will throw an exception. We therefore set here the use_model
901  // flag to false in case that there are no spectral nodes
902  if (nodes->nodes() == 0) {
903  use_model = false;
904  }
905 
906  } // endif: spatial model was a diffuse cube
907 
908  // Compute flux within [emin, emax] in model from spectral
909  // component (units: ph/cm2/s)
910  double flux = (use_model) ? spectral->flux(emin, emax) : 0.0;
911 
912  // Derive expecting counting rate within simulation surface
913  // (units: ph/s)
914  double rate = flux * area * norm;
915 
916  // Debug option: dump rate
917  #if defined(G_DUMP_MC)
918  std::cout << "GModelSky::mc(\"" << name() << "\": ";
919  std::cout << "flux=" << flux << " ph/cm2/s, ";
920  std::cout << "rate=" << rate << " ph/s, ";
921  std::cout << "norm=" << norm << ")" << std::endl;
922  #endif
923 
924  // Continue only if rate is positive
925  if (rate > 0.0) {
926 
927  // Get photon arrival times from temporal model
928  GTimes times = m_temporal->mc(rate, tmin, tmax, ran);
929 
930  // Debug option: dump number of times
931  #if defined(G_DUMP_MC_DETAIL)
932  std::cout << " Times=" << times.size() << std::endl;
933  #endif
934 
935  // Reserve space for photons
936  if (times.size() > 0) {
937  photons.reserve(times.size());
938  }
939 
940  // Loop over photons
941  for (int i = 0; i < times.size(); ++i) {
942 
943  // Debug option: dump photon index
944  #if defined(G_DUMP_MC_DETAIL)
945  std::cout << " Photon=" << i << std::endl;
946  #endif
947 
948  // Allocate photon
949  GPhoton photon;
950 
951  // Set photon arrival time
952  photon.time(times[i]);
953 
954  // Debug option: dump time
955  #if defined(G_DUMP_MC_DETAIL)
956  std::cout << " Time=" << times[i] << std::endl;
957  #endif
958 
959  // Set photon energy. If an invalid_return_value exception
960  // occurs the energy returned by the spectral Monte Carlo
961  // method is invalid and the photon is skipped.
962  try {
963  photon.energy(spectral->mc(emin, emax, photon.time(),
964  ran));
965  }
967  continue;
968  }
969 
970  // Debug option: dump energy
971  #if defined(G_DUMP_MC_DETAIL)
972  std::cout << " Energy=" << photon.energy() << std::endl;
973  #endif
974 
975  // Set incident photon direction. If an invalid_return_value
976  // exception occurs the sky direction returned by the
977  // spatial Monte Carlo method is invalid and the photon is
978  // skipped.
979  try {
980  photon.dir(m_spatial->mc(photon.energy(), photon.time(),
981  ran));
982  }
984  continue;
985  }
986 
987  // Debug option: dump direction
988  #if defined(G_DUMP_MC_DETAIL)
989  std::cout << " Direction=" << photon.dir() << std::endl;
990  #endif
991 
992  // Append photon
993  if (dir.dist_deg(photon.dir()) <= radius) {
994  photons.append(photon);
995  }
996 
997  } // endfor: looped over photons
998 
999  } // endif: rate was positive
1000 
1001  // Free spectral model if required
1002  if (free_spectral) delete spectral;
1003 
1004  } // endif: model was used
1005 
1006  } // endif: model was valid
1007 
1008  // Return photon list
1009  return photons;
1010 }
1011 
1012 
1013 /***********************************************************************//**
1014  * @brief Print model information
1015  *
1016  * @param[in] chatter Chattiness.
1017  * @return String containing model information.
1018  ***************************************************************************/
1019 std::string GModelSky::print(const GChatter& chatter) const
1020 {
1021  // Initialise result string
1022  std::string result;
1023 
1024  // Continue only if chatter is not silent
1025  if (chatter != SILENT) {
1026 
1027  // Append header
1028  result.append("=== GModelSky ===");
1029 
1030  // Determine number of parameters per type
1031  int n_spatial = (m_spatial != NULL) ? m_spatial->size() : 0;
1032  int n_spectral = (m_spectral != NULL) ? m_spectral->size() : 0;
1033  int n_temporal = (m_temporal != NULL) ? m_temporal->size() : 0;
1034 
1035  // Append attributes
1036  result.append("\n"+print_attributes());
1037 
1038  // Append model type
1039  result.append("\n"+gammalib::parformat("Model type")+type());
1040 
1041  // Append model components
1042  result.append("\n"+gammalib::parformat("Model components"));
1043  if (n_spatial > 0) {
1044  result.append("\""+spatial()->type()+"\"");
1045  if (n_spectral > 0 || n_temporal > 0) {
1046  result.append(" * ");
1047  }
1048  }
1049  if (n_spectral > 0) {
1050  result.append("\""+spectral()->type()+"\"");
1051  if (n_temporal > 0) {
1052  result.append(" * ");
1053  }
1054  }
1055  if (n_temporal > 0) {
1056  result.append("\""+temporal()->type()+"\"");
1057  }
1058 
1059  // Append parameters
1060  result.append("\n"+gammalib::parformat("Number of parameters"));
1061  result.append(gammalib::str(size()));
1062  result.append("\n"+gammalib::parformat("Number of spatial par's"));
1063  result.append(gammalib::str(n_spatial));
1064  for (int i = 0; i < n_spatial; ++i) {
1065  result.append("\n"+(*spatial())[i].print());
1066  }
1067  result.append("\n"+gammalib::parformat("Number of spectral par's"));
1068  result.append(gammalib::str(n_spectral));
1069  for (int i = 0; i < n_spectral; ++i) {
1070  result.append("\n"+(*spectral())[i].print());
1071  }
1072  result.append("\n"+gammalib::parformat("Number of temporal par's"));
1073  result.append(gammalib::str(n_temporal));
1074  for (int i = 0; i < n_temporal; ++i) {
1075  result.append("\n"+(*temporal())[i].print());
1076  }
1077  result.append("\n"+gammalib::parformat("Number of scale par's"));
1078  result.append(gammalib::str(scales()));
1079  for (int i = 0; i < scales(); ++i) {
1080  result.append("\n"+scale(i).print());
1081  }
1082 
1083  } // endif: chatter was not silent
1084 
1085  // Return result
1086  return result;
1087 }
1088 
1089 
1090 /*==========================================================================
1091  = =
1092  = Private methods =
1093  = =
1094  ==========================================================================*/
1095 
1096 /***********************************************************************//**
1097  * @brief Initialise class members
1098  ***************************************************************************/
1100 {
1101  // Initialise members
1102  m_type.clear();
1103  m_spatial = NULL;
1104  m_spectral = NULL;
1105  m_temporal = NULL;
1106 
1107  // Return
1108  return;
1109 }
1110 
1111 
1112 /***********************************************************************//**
1113  * @brief Copy class members
1114  *
1115  * @param[in] model Sky model.
1116  ***************************************************************************/
1118 {
1119  // Copy attributes
1120  m_type = model.m_type;
1121 
1122  // Clone model components
1123  m_spatial = (model.m_spatial != NULL) ? model.m_spatial->clone() : NULL;
1124  m_spectral = (model.m_spectral != NULL) ? model.m_spectral->clone() : NULL;
1125  m_temporal = (model.m_temporal != NULL) ? model.m_temporal->clone() : NULL;
1126 
1127  // Set parameter pointers
1128  set_pointers();
1129 
1130  // Return
1131  return;
1132 }
1133 
1134 
1135 /***********************************************************************//**
1136  * @brief Delete class members
1137  ***************************************************************************/
1139 {
1140  // Free memory
1141  if (m_spatial != NULL) delete m_spatial;
1142  if (m_spectral != NULL) delete m_spectral;
1143  if (m_temporal != NULL) delete m_temporal;
1144 
1145  // Signal free pointers
1146  m_spatial = NULL;
1147  m_spectral = NULL;
1148  m_temporal = NULL;
1149 
1150  // Return
1151  return;
1152 }
1153 
1154 
1155 /***********************************************************************//**
1156  * @brief Set parameter pointers
1157  *
1158  * Gathers all parameter pointers from the model into a flat array of model
1159  * pointers.
1160  ***************************************************************************/
1162 {
1163  // Clear parameters
1164  m_pars.clear();
1165 
1166  // Determine the number of parameters
1167  int n_spatial = (spatial() != NULL) ? spatial()->size() : 0;
1168  int n_spectral = (spectral() != NULL) ? spectral()->size() : 0;
1169  int n_temporal = (temporal() != NULL) ? temporal()->size() : 0;
1170  int n_pars = n_spatial + n_spectral + n_temporal;
1171  int n_scales = m_scales.size();
1172 
1173  // Continue only if there are parameters
1174  if (n_pars > 0) {
1175 
1176  // Gather spatial parameter pointers
1177  for (int i = 0; i < n_spatial; ++i) {
1178  m_pars.push_back(&((*spatial())[i]));
1179  }
1180 
1181  // Gather spectral parameters
1182  for (int i = 0; i < n_spectral; ++i) {
1183  m_pars.push_back(&((*spectral())[i]));
1184  }
1185 
1186  // Gather temporal parameters
1187  for (int i = 0; i < n_temporal; ++i) {
1188  m_pars.push_back(&((*temporal())[i]));
1189  }
1190 
1191  }
1192 
1193  // Append scales if they exist
1194  if (n_scales > 0) {
1195  for (int i = 0; i < n_scales; ++i) {
1196  m_pars.push_back(&m_scales[i]);
1197  }
1198  }
1199 
1200  // Return
1201  return;
1202 }
1203 
1204 
1205 /***********************************************************************//**
1206  * @brief Set model type based on spatial model component
1207  *
1208  * Set the model type depending on the class type of the spatial model
1209  * component.
1210  *
1211  * @todo A method could be implemented in the GModelSpatial class that
1212  * determines the model type. This is however not very critical.
1213  ***************************************************************************/
1215 {
1216  // Clear model type
1217  m_type.clear();
1218 
1219  // Continue only if we have a spatial model component
1220  if (m_spatial != NULL) {
1221 
1222  // Is spatial model a point source?
1223  if (dynamic_cast<const GModelSpatialPointSource*>(m_spatial) != NULL) {
1224  m_type = "PointSource";
1225  }
1226 
1227  // ... otherwise, is spatial model a radial source?
1228  else if (dynamic_cast<const GModelSpatialRadial*>(m_spatial) != NULL) {
1229  m_type = "ExtendedSource";
1230  }
1231 
1232  // ... otherwise, is spatial model an elliptical source?
1233  else if (dynamic_cast<const GModelSpatialElliptical*>(m_spatial) != NULL) {
1234  m_type = "ExtendedSource";
1235  }
1236 
1237  // ... otherwise, is spatial model a composite source?
1238  else if (dynamic_cast<const GModelSpatialComposite*>(m_spatial) != NULL) {
1239  m_type = "CompositeSource";
1240  }
1241 
1242  // ... otherwise we have a diffuse model
1243  else {
1244  m_type = "DiffuseSource";
1245  }
1246 
1247  } // endif: there was a spatial model component
1248 
1249  // Return
1250  return;
1251 }
1252 
1253 
1254 /***********************************************************************//**
1255  * @brief Return pointer to spatial model from XML element
1256  *
1257  * @param[in] spatial XML element.
1258  * @return Pointer to spatial model.
1259  *
1260  * Returns pointer to spatial model that is defined in an XML element.
1261  ***************************************************************************/
1263 {
1264  // Get spatial model
1265  GModelSpatialRegistry registry;
1266  GModelSpatial* ptr = registry.alloc(spatial);
1267 
1268  // Return pointer
1269  return ptr;
1270 }
1271 
1272 
1273 /***********************************************************************//**
1274  * @brief Return pointer to spectral model from XML element
1275  *
1276  * @param[in] spectral XML element.
1277  * @return Pointer to spectral model.
1278  *
1279  * Returns pointer to spectral model that is defined in an XML element.
1280  ***************************************************************************/
1282 {
1283  // Get spectral model
1284  GModelSpectralRegistry registry;
1285  GModelSpectral* ptr = registry.alloc(spectral);
1286 
1287  // Return pointer
1288  return ptr;
1289 }
1290 
1291 
1292 /***********************************************************************//**
1293  * @brief Return pointer to temporal model from XML element
1294  *
1295  * @param[in] temporal XML element.
1296  * @return Pointer to temporal model.
1297  *
1298  * Returns pointer to temporal model that is defined in an XML element.
1299  ***************************************************************************/
1301 {
1302  // Get temporal model
1303  GModelTemporalRegistry registry;
1304  GModelTemporal* ptr = registry.alloc(temporal);
1305 
1306  // Return pointer
1307  return ptr;
1308 }
1309 
1310 
1311 /***********************************************************************//**
1312  * @brief Verifies if model has all components
1313  *
1314  * @return True is model is valid, false otherwise.
1315  ***************************************************************************/
1316 bool GModelSky::valid_model(void) const
1317 {
1318  // Set result
1319  bool result = ((m_spatial != NULL) &&
1320  (m_spectral != NULL) &&
1321  (m_temporal != NULL));
1322 
1323  // Return result
1324  return result;
1325 }
virtual void write(GXmlElement &xml) const =0
Constant temporal model class interface definition.
Abstract model class.
Definition: GModel.hpp:97
double dist_deg(const GSkyDir &dir) const
Compute angular distance between sky directions in degrees.
Definition: GSkyDir.hpp:280
void free_members(void)
Delete class members.
Definition: GModel.cpp:656
int scales(void) const
Return number of scale parameters in model.
Definition: GModel.hpp:231
void set_type(void)
Set model type based on spatial model component.
Definition: GModelSky.cpp:1214
virtual GSkyDir mc(const GEnergy &energy, const GTime &time, GRan &ran) const =0
GModelTemporal * m_temporal
Temporal model.
Definition: GModelSky.hpp:190
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:821
Point source spatial model class interface definition.
int size(void) const
Return number of times.
Definition: GTimes.hpp:100
void set_mc_cone(const GSkyDir &centre, const double &radius) const
Set Monte Carlo simulation cone.
virtual double eval(const GPhoton &photon, const bool &gradients=false) const =0
std::string print_attributes(void) const
Print model attributes.
Definition: GModel.cpp:748
Abstract spectral model base class.
void copy_members(const GModelSky &model)
Copy class members.
Definition: GModelSky.cpp:1117
Interface definition for the model registry class.
GModelSpectral * spectral(void) const
Return spectral model component.
Definition: GModelSky.hpp:262
GModelTemporal * xml_temporal(const GXmlElement &temporal) const
Return pointer to temporal model from XML element.
Definition: GModelSky.cpp:1300
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:845
int size(void) const
Return number of parameters.
virtual void clear(void)
Clear sky model.
Definition: GModelSky.cpp:378
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:1262
GModelTemporal * temporal(void) const
Return temporal model component.
Definition: GModelSky.hpp:278
Abstract radial spatial model base class interface definition.
Abstract interface for the event classes.
Definition: GEvent.hpp:71
void set_pointers(void)
Set parameter pointers.
Definition: GModelSky.cpp:1161
XML element node class.
Definition: GXmlElement.hpp:47
Spectral model registry class definition.
virtual ~GModelSky(void)
Destructor.
Definition: GModelSky.cpp:321
virtual void write(GXmlElement &xml) const =0
Random number generator class.
Definition: GRan.hpp:44
Photon container class.
Definition: GPhotons.hpp:45
bool valid_model(void) const
Verifies if model has all components.
Definition: GModelSky.cpp:1316
virtual int elements(void) const
Return number of GXMLElement children of node.
Definition: GXmlNode.cpp:580
Time class.
Definition: GTime.hpp:54
Gammalib tools definition.
virtual GModelTemporalConst * clone(void) const
Clone constant temporal model.
Interface definition for the spatial model registry class.
Source class definition.
std::string m_type
Model type.
Definition: GModelSky.hpp:187
GModelSpectral * m_spectral
Spectral model.
Definition: GModelSky.hpp:189
std::vector< GModelPar > m_scales
Model instrument scale factors.
Definition: GModel.hpp:167
Spectral nodes model class.
Time container class.
Definition: GTimes.hpp:45
bool is_notanumber(const double &x)
Signal if argument is not a number.
Definition: GTools.hpp:184
std::vector< GModelPar * > m_pars
Pointers to all model parameters.
Definition: GModel.hpp:169
virtual void write(GXmlElement &xml) const
Write model into XML element.
Definition: GModelSky.cpp:756
int size(void) const
Return number of parameters in model.
Definition: GModel.hpp:217
Class that handles photons.
Definition: GPhoton.hpp:47
bool is_infinite(const double &x)
Signal if argument is infinite.
Definition: GTools.hpp:167
void init_members(void)
Initialise class members.
Definition: GModel.cpp:612
const GModelSpectralNodes & spectrum(void) const
Get map cube spectrum.
const GXmlAttribute * attribute(const int &index) const
Return attribute.
const std::string & name(void) const
Return parameter name.
Definition: GModel.hpp:245
GVector gradients(const GPhoton &photon)
Return parameter gradients of sky model for a given photon.
Definition: GModelSky.cpp:528
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:611
GModelTemporal * alloc(const GXmlElement &xml) const
Allocate temporal model that is found in XML element.
Spatial model registry class definition.
virtual void read(const GXmlElement &xml)
Read sky model from XML element.
Definition: GModelSky.cpp:687
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
Spatial composite model class interface definition.
Abstract observation base class.
virtual GModelSky * clone(void) const
Clone sky model.
Definition: GModelSky.cpp:398
const GModelSky g_extendedsource_seed("ExtendedSource")
Definition: GModelSky.cpp:54
GModelPar & scale(const int &index)
Returns reference to scale parameter by index.
Definition: GModel.cpp:363
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:704
virtual GModelSpectral * clone(void) const =0
Clones object.
virtual GModelSky & operator=(const GModelSky &model)
Assignment operator.
Definition: GModelSky.cpp:343
Abstract elliptical spatial model base class interface definition.
Abstract response base class definition.
double intensity(const int &index) const
Return node intensity.
double value(const GPhoton &photon)
Return value of sky model for a given photon.
Definition: GModelSky.cpp:502
const GModelSky g_diffusesource_seed("DiffuseSource")
Definition: GModelSky.cpp:55
virtual GXmlElement * element(const int &index)
Return pointer to GXMLElement child.
Definition: GXmlNode.cpp:634
void reserve(const int &num)
Reserve memory for photons in container.
Definition: GPhotons.cpp:327
Temporal model registry class definition.
Sky model class interface definition.
Sky model class.
Definition: GModelSky.hpp:120
Spatial map cube model class interface definition.
int nodes(void) const
Return number of nodes.
void init_members(void)
Initialise class members.
Definition: GModelSky.cpp:1099
virtual double mc_norm(const GSkyDir &dir, const double &radius) const =0
void append(const GPhoton &photon)
Append photon to container.
Definition: GPhotons.cpp:215
virtual std::string type(void) const
Return sky model type.
Definition: GModelSky.hpp:215
GModelSpatial * m_spatial
Spatial model.
Definition: GModelSky.hpp:188
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:246
virtual GModelTemporal * clone(void) const =0
Clones object.
void read_attributes(const GXmlElement &xml)
Read model attributes.
Definition: GModel.cpp:668
Abstract spatial model base class.
virtual void write(GXmlElement &xml) const =0
GModelSky(void)
Void constructor.
Definition: GModelSky.cpp:90
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 GXmlNode * append(const GXmlNode &node)
Append XML child node.
Definition: GXmlNode.cpp:284
std::string parformat(const std::string &s, const int &indent=0)
Convert string in parameter format.
Definition: GTools.cpp:1022
virtual std::string print(const GChatter &chatter=NORMAL) const
Print model information.
Definition: GModelSky.cpp:1019
Integration class interface definition.
Sky direction class.
Definition: GSkyDir.hpp:62
const GModelSky g_compositesource_seed("CompositeSource")
Definition: GModelSky.cpp:56
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:564
const GModelSky g_pointsource_seed("PointSource")
Definition: GModelSky.cpp:53
Constant temporal model class.
void free_members(void)
Delete class members.
Definition: GModelSky.cpp:1138
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:1281
std::string str(const unsigned short int &value)
Convert unsigned short integer value into string.
Definition: GTools.cpp:413