GammaLib  1.7.0.dev
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
GCTAModelBackground.cpp
Go to the documentation of this file.
1 /***************************************************************************
2  * GCTAModelBackground.cpp - Background model class *
3  * ----------------------------------------------------------------------- *
4  * copyright (C) 2018-2020 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 GCTAModelBackground.cpp
23  * @brief Background model class implementation
24  * @author Juergen Knoedlseder
25  */
26 
27 /* __ Includes ___________________________________________________________ */
28 #ifdef HAVE_CONFIG_H
29 #include <config.h>
30 #endif
31 #include "GException.hpp"
32 #include "GTools.hpp"
33 #include "GMath.hpp"
34 #include "GModelRegistry.hpp"
37 #include "GModelTemporalConst.hpp"
38 #include "GIntegral.hpp"
40 #include "GCTAModelBackground.hpp"
41 #include "GCTAObservation.hpp"
42 #include "GCTAPointing.hpp"
43 #include "GCTAInstDir.hpp"
44 #include "GCTARoi.hpp"
45 #include "GCTAException.hpp"
46 #include "GCTASupport.hpp"
47 
48 /* __ Constants __________________________________________________________ */
49 
50 /* __ Globals ____________________________________________________________ */
52 const GModelRegistry g_cta_background_registry(&g_cta_background_seed);
53 
54 /* __ Method name definitions ____________________________________________ */
55 #define G_EVAL "GCTAModelBackground::eval(GEvent&, GObservation&, bool&)"
56 #define G_NPRED "GCTAModelBackground::npred(GEnergy&, GTime&, GObservation&)"
57 #define G_MC "GCTAModelBackground::mc(GObservation&, GRan&)"
58 #define G_XML_SPATIAL "GCTAModelBackground::xml_spatial(GXmlElement&)"
59 #define G_XML_SPECTRAL "GCTAModelBackground::xml_spectral(GXmlElement&)"
60 #define G_XML_TEMPORAL "GCTAModelBackground::xml_temporal(GXmlElement&)"
61 
62 /* __ Macros _____________________________________________________________ */
63 
64 /* __ Coding definitions _________________________________________________ */
65 
66 /* __ Debug definitions __________________________________________________ */
67 //#define G_DUMP_MC //!< Dump MC information
68 
69 
70 /*==========================================================================
71  = =
72  = Constructors/destructors =
73  = =
74  ==========================================================================*/
75 
76 /***********************************************************************//**
77  * @brief Void constructor
78  *
79  * Constructs an empty background model.
80  ***************************************************************************/
82 {
83  // Initialise members
84  init_members();
85 
86  // Return
87  return;
88 }
89 
90 
91 /***********************************************************************//**
92  * @brief XML constructor
93  *
94  * @param[in] xml XML element.
95  *
96  * Constructs a background model from the information that is found in a XML
97  * element. Please refer to the read() method to learn about the expected
98  * structure of the XML element.
99  ***************************************************************************/
101  GModelData(xml)
102 {
103  // Initialise members
104  init_members();
105 
106  // Read XML
107  read(xml);
108 
109  // Set parameter pointers
110  set_pointers();
111 
112  // Return
113  return;
114 }
115 
116 
117 /***********************************************************************//**
118  * @brief Construct from spatial and spectral components
119  *
120  * @param[in] spatial Spatial model component.
121  * @param[in] spectral Spectral model component.
122  *
123  * Constructs a background model from a @p spatial and a @p spectral model
124  * component. The temporal component is assumed to be constant.
125  ***************************************************************************/
127  const GModelSpectral& spectral) :
128  GModelData()
129 {
130  // Initialise members
131  init_members();
132 
133  // Allocate temporal constant model
135 
136  // Clone model components
137  m_spatial = spatial.clone();
138  m_spectral = spectral.clone();
139  m_temporal = temporal.clone();
140 
141  // Set parameter pointers
142  set_pointers();
143 
144  // Return
145  return;
146 }
147 
148 
149 /***********************************************************************//**
150  * @brief Construct from model components
151  *
152  * @param[in] spatial Spatial model component.
153  * @param[in] spectral Spectral model component.
154  * @param[in] temporal Temporal model component.
155  *
156  * Constructs a background model from a @p spatial, a @p spectral and a
157  * @p temporal model component.
158  ***************************************************************************/
160  const GModelSpectral& spectral,
161  const GModelTemporal& temporal) :
162  GModelData()
163 {
164  // Initialise members
165  init_members();
166 
167  // Clone model components
168  m_spatial = spatial.clone();
169  m_spectral = spectral.clone();
170  m_temporal = temporal.clone();
171 
172  // Set parameter pointers
173  set_pointers();
174 
175  // Return
176  return;
177 }
178 
179 
180 /***********************************************************************//**
181  * @brief Copy constructor
182  *
183  * @param[in] model Background model.
184  *
185  * Constructs a background model by copying information from an existing
186  * model. Note that the copy is a deep copy, so the original object can be
187  * destroyed after the copy without any loss of information.
188  ***************************************************************************/
190  GModelData(model)
191 {
192  // Initialise private members for clean destruction
193  init_members();
194 
195  // Copy members
196  copy_members(model);
197 
198  // Set parameter pointers
199  set_pointers();
200 
201  // Return
202  return;
203 }
204 
205 
206 /***********************************************************************//**
207  * @brief Destructor
208  *
209  * Destroys a background model.
210  ***************************************************************************/
212 {
213  // Free members
214  free_members();
215 
216  // Return
217  return;
218 }
219 
220 
221 /*==========================================================================
222  = =
223  = Operators =
224  = =
225  ==========================================================================*/
226 
227 /***********************************************************************//**
228  * @brief Assignment operator
229  *
230  * @param[in] model Background model.
231  *
232  * Assigns the information from a background model to instance. Note that a
233  * deep copy of the information is performed, so the original instance can be
234  * destroyed after the assignment without any loss of information.
235  ***************************************************************************/
237 {
238  // Execute only if object is not identical
239  if (this != &model) {
240 
241  // Copy base class members
242  this->GModelData::operator=(model);
243 
244  // Free members
245  free_members();
246 
247  // Initialise members
248  init_members();
249 
250  // Copy members (this method also sets the parameter pointers)
251  copy_members(model);
252 
253  } // endif: object was not identical
254 
255  // Return
256  return *this;
257 }
258 
259 
260 /*==========================================================================
261  = =
262  = Public methods =
263  = =
264  ==========================================================================*/
265 
266 /***********************************************************************//**
267  * @brief Clear instance
268  *
269  * Resets the instance to a clean initial state. All information that resided
270  * in the object will be lost.
271  ***************************************************************************/
273 {
274  // Free class members (base and derived classes, derived class first)
275  free_members();
276  this->GModelData::free_members();
277  this->GModel::free_members();
278 
279  // Initialise members
280  this->GModel::init_members();
281  this->GModelData::init_members();
282  init_members();
283 
284  // Return
285  return;
286 }
287 
288 
289 /***********************************************************************//**
290  * @brief Clone instance
291  *
292  * @return Pointer to deep copy of background model.
293  ***************************************************************************/
295 {
296  return new GCTAModelBackground(*this);
297 }
298 
299 
300 /***********************************************************************//**
301  * @brief Set spatial model component
302  *
303  * @param[in] spatial Pointer to spatial model component.
304  *
305  * Sets the spatial model component of the model.
306  ***************************************************************************/
308 {
309  // Free spatial model component only if it differs from current
310  // component. This prevents unintentional deallocation of the argument
311  if ((m_spatial != NULL) && (m_spatial != spatial)) {
312  delete m_spatial;
313  }
314 
315  // Clone spatial model component if it exists, otherwise set pointer
316  // to NULL
317  m_spatial = (spatial != NULL) ? spatial->clone() : NULL;
318 
319  // Set parameter pointers
320  set_pointers();
321 
322  // Return
323  return;
324 }
325 
326 
327 /***********************************************************************//**
328  * @brief Set spectral model component
329  *
330  * @param[in] spectral Pointer to spectral model component.
331  *
332  * Sets the spectral model component of the model.
333  ***************************************************************************/
335 {
336  // Free spectral model component only if it differs from current
337  // component. This prevents unintential deallocation of the argument
338  if ((m_spectral != NULL) && (m_spectral != spectral)) {
339  delete m_spectral;
340  }
341 
342  // Clone spectral model component if it exists, otherwise set pointer
343  // to NULL
344  m_spectral = (spectral != NULL) ? spectral->clone() : NULL;
345 
346  // Set parameter pointers
347  set_pointers();
348 
349  // Return
350  return;
351 }
352 
353 
354 /***********************************************************************//**
355  * @brief Set temporal model component
356  *
357  * @param[in] temporal Pointer to temporal model component.
358  *
359  * Sets the temporal model component of the model.
360  ***************************************************************************/
362 {
363  // Free temporal model component only if it differs from current
364  // component. This prevents unintential deallocation of the argument
365  if ((m_temporal != NULL) && (m_temporal != temporal)) {
366  delete m_temporal;
367  }
368 
369  // Clone temporal model component if it exists, otherwise set pointer
370  // to NULL
371  m_temporal = (temporal != NULL) ? temporal->clone() : NULL;
372 
373  // Set parameter pointers
374  set_pointers();
375 
376  // Return
377  return;
378 }
379 
380 
381 /***********************************************************************//**
382  * @brief Return background rate in units of events MeV\f$^{-1}\f$
383  * s\f$^{-1}\f$ sr\f$^{-1}\f$
384  *
385  * @param[in] event Observed event.
386  * @param[in] obs Observation.
387  * @param[in] gradients Compute gradients?
388  * @return Background rate (events MeV\f$^{-1}\f$ s\f$^{-1}\f$ sr\f$^{-1}\f$).
389  *
390  * Evaluates the background model which is a factorization of a spatial,
391  * spectral and temporal model component. The method returns a real rate,
392  * defined as the number of counts per MeV, steradian and ontime.
393  *
394  * If the @p gradients flag is true the method will also set the parameter
395  * gradients of the model parameters.
396  *
397  * @todo Add bookkeeping of last value and evaluate only if argument
398  * changed
399  * @todo Verify that CTA instrument direction pointer is valid, or better,
400  * add an offset method to GCTAPointing. Ideally, we should precompute
401  * all offset angles (for an event cube this may only be done easily
402  * if the pointing has been fixed; otherwise we need a structure
403  * similar to the Fermi/LAT livetime cube that provides the effective
404  * sky exposure as function of offset angle).
405  ***************************************************************************/
406 double GCTAModelBackground::eval(const GEvent& event,
407  const GObservation& obs,
408  const bool& gradients) const
409 {
410  // Get reference on CTA instrument direction from event
411  const GCTAInstDir& dir = gammalib::cta_dir(G_EVAL, event);
412 
413  // Evaluate function and gradients
414  double spat = (spatial() != NULL)
415  ? spatial()->eval(dir, event.energy(), event.time(), gradients) : 1.0;
416  double spec = (spectral() != NULL)
417  ? spectral()->eval(event.energy(), event.time(), gradients)
418  : 1.0;
419  double temp = (temporal() != NULL)
420  ? temporal()->eval(event.time(), gradients) : 1.0;
421 
422  // Compute value
423  double value = spat * spec * temp;
424 
425  // Optionally compute partial derivatives
426  if (gradients) {
427 
428  // Multiply factors to radial gradients
429  if (spatial() != NULL) {
430  double fact = spec * temp;
431  if (fact != 1.0) {
432  for (int i = 0; i < spatial()->size(); ++i)
433  (*spatial())[i].factor_gradient((*spatial())[i].factor_gradient() * fact);
434  }
435  }
436 
437  // Multiply factors to spectral gradients
438  if (spectral() != NULL) {
439  double fact = spat * temp;
440  if (fact != 1.0) {
441  for (int i = 0; i < spectral()->size(); ++i)
442  (*spectral())[i].factor_gradient((*spectral())[i].factor_gradient() * fact);
443  }
444  }
445 
446  // Multiply factors to temporal gradients
447  if (temporal() != NULL) {
448  double fact = spat * spec;
449  if (fact != 1.0) {
450  for (int i = 0; i < temporal()->size(); ++i)
451  (*temporal())[i].factor_gradient((*temporal())[i].factor_gradient() * fact);
452  }
453  }
454 
455  } // endif: computed partial derivatives
456 
457  // Return
458  return value;
459 }
460 
461 
462 /***********************************************************************//**
463  * @brief Return spatially integrated background rate in units of
464  * events MeV\f$^{-1}\f$ s\f$^{-1}\f$
465  *
466  * @param[in] energy Measured event energy.
467  * @param[in] time Measured event time.
468  * @param[in] obs Observation.
469  * @return Spatially integrated background rate
470  * (events MeV\f$^{-1}\f$ s\f$^{-1}\f$)
471  *
472  * Spatially integrates the background model for a given measured event
473  * energy and event time. The method returns a real rate, defined as the
474  * number of counts per MeV and ontime.
475  ***************************************************************************/
476 double GCTAModelBackground::npred(const GEnergy& energy,
477  const GTime& time,
478  const GObservation& obs) const
479 {
480  // Get spatially integrated model component
481  double npred = spatial()->npred(energy, time, obs);
482 
483  // Multiply in spectral and temporal components
484  npred *= spectral()->eval(energy, time);
485  npred *= temporal()->eval(time);
486 
487  // Return Npred
488  return npred;
489 }
490 
491 
492 /***********************************************************************//**
493  * @brief Return simulated list of events
494  *
495  * @param[in] obs Observation.
496  * @param[in] ran Random number generator.
497  * @return Pointer to list of simulated events (needs to be de-allocated by
498  * client)
499  *
500  * Draws a sample of events from the background model using a Monte
501  * Carlo simulation. The region of interest, the energy boundaries and the
502  * good time interval for the sampling will be extracted from the observation
503  * argument that is passed to the method. The method also requires a random
504  * number generator of type GRan which is passed by reference, hence the
505  * state of the random number generator will be changed by the method.
506  *
507  * For each event in the returned event list, the sky direction, the nominal
508  * coordinates (DETX and DETY), the energy and the time will be set.
509  *
510  * @todo Handle special case of cube spatial model
511  ***************************************************************************/
513  GRan& ran) const
514 {
515  // Initialise new event list
516  GCTAEventList* list = new GCTAEventList;
517 
518  // Continue only if model is valid)
519  if (valid_model()) {
520 
521  // Get reference on CTA pointing, background response and event list
522  // from observation
523  const GCTAObservation& cta = gammalib::cta_obs(G_MC, obs);
524  const GCTAPointing& pnt = gammalib::cta_pnt(G_MC, obs);
525  const GCTAEventList& events = gammalib::cta_event_list(G_MC, obs);
526 
527  // Get simulation region
528  const GCTARoi& roi = events.roi();
529  const GEbounds& ebounds = events.ebounds();
530  const GGti& gti = events.gti();
531 
532  // Set simulation region for result event list
533  list->roi(roi);
534  list->ebounds(ebounds);
535  list->gti(gti);
536 
537  // Set spectral model
539 
540  // TODO: Handle special case of cube spatial model. This will
541  // replace the spectral model by a spectral nodes model.
542 
543  // Get solid angle of spatial model. This only works for an energy
544  // independent spatial model!!!
545  double solidangle = m_spatial->npred(GEnergy(), GTime(), obs);
546 
547  // Loop over all energy boundaries
548  for (int ieng = 0; ieng < ebounds.size(); ++ieng) {
549 
550  // Compute the background rate in model within the energy
551  // boundaries from spectral component (units: cts/s/sr).
552  double flux = spectral->flux(ebounds.emin(ieng), ebounds.emax(ieng));
553 
554  // Derive expecting rate (units: cts/s). Note that the time here
555  // is good time. Deadtime correction will be done later.
556  double rate = flux * solidangle;
557 
558  // Debug option: dump rate
559  #if defined(G_DUMP_MC)
560  std::cout << "GCTAModelBackground::mc(\"" << name() << "\": ";
561  std::cout << "flux=" << flux << " cts/s/sr, ";
562  std::cout << "solidangle=" << solidangle << " sr, ";
563  std::cout << "rate=" << rate << " cts/s)" << std::endl;
564  #endif
565 
566  // If the rate is not positive then skip this energy bins
567  if (rate <= 0.0) {
568  continue;
569  }
570 
571  // Loop over all good time intervals
572  for (int itime = 0; itime < gti.size(); ++itime) {
573 
574  // Get Monte Carlo event arrival times from temporal model
575  GTimes times = m_temporal->mc(rate,
576  gti.tstart(itime),
577  gti.tstop(itime),
578  ran);
579 
580  // Get number of events
581  int n_events = times.size();
582 
583  // Reserve space for events
584  if (n_events > 0) {
585  list->reserve(n_events);
586  }
587 
588  // Debug option: provide number of times and initialize
589  // statisics
590  #if defined(G_DUMP_MC)
591  std::cout << " Interval " << itime;
592  std::cout << " events=" << n_events << std::endl;
593  int n_trial_outside_roi = 0;
594  #endif
595 
596  // Loop over events
597  for (int i = 0; i < n_events; ++i) {
598 
599  // Get Monte Carlo event energy from spectral model
600  GEnergy energy = spectral->mc(ebounds.emin(ieng),
601  ebounds.emax(ieng),
602  times[i],
603  ran);
604 
605  // Get an instrument direction within the RoI. This is
606  // potentially tried 100 times so that if we really can't
607  // a valid instrument direction the code is not locked up
608  for (int k = 0; k < 100; ++k) {
609 
610  // Get Monte Carlo event direction from spatial model.
611  GCTAInstDir dir = m_spatial->mc(energy, times[i], cta, ran);
612 
613  // Allocate event
614  GCTAEventAtom event;
615 
616  // Set event attributes
617  event.dir(dir);
618  event.energy(energy);
619  event.time(times[i]);
620 
621  // Append event to list if it falls in RoI and break
622  // the look
623  if (events.roi().contains(event)) {
624  list->append(event);
625  break;
626  }
627 
628  // ... otherwise optionally bookkeep the trial
629  #if defined(G_DUMP_MC)
630  else {
631  n_trial_outside_roi++;
632  }
633  #endif
634 
635  } // endfor: trial look for instrument direction within RoI
636 
637  } // endfor: looped over all events
638 
639  // Debug option: provide statisics
640  #if defined(G_DUMP_MC)
641  std::cout << " Trials outside RoI=";
642  std::cout << n_trial_outside_roi << std::endl;
643  #endif
644 
645  } // endfor: looped over all GTIs
646 
647  } // endfor: looped over all energy boundaries
648 
649  } // endif: model was valid
650 
651  // Return
652  return list;
653 }
654 
655 
656 /***********************************************************************//**
657  * @brief Read model from XML element
658  *
659  * @param[in] xml XML element.
660  *
661  * Reads the sky model from an XML element. The XML element is expected to
662  * respect the following format:
663  *
664  * <source name=".." type=".." instrument=".." id="..">
665  * <spectrum type="..">
666  * ..
667  * </spectrum>
668  * <spatialModel type="..">
669  * ..
670  * </spatialModel>
671  * <temporal type="..">
672  * ..
673  * </temporal>
674  * </source>
675  *
676  * The temporal element is optional. In no temporal element is specified a
677  * constant component with unity normalization will be assumed.
678  ***************************************************************************/
680 {
681  // Clear model
682  clear();
683 
684  // Get pointers on spatial and spectral model components
685  const GXmlElement* spat = xml.element("spatialModel", 0);
686  const GXmlElement* spec = xml.element("spectrum", 0);
687 
688  // Set spatial and spectral model components
689  m_spatial = xml_spatial(*spat);
690  m_spectral = xml_spectral(*spec);
691 
692  // Handle optional temporal model
693  if (xml.elements("temporal") > 0) {
694  const GXmlElement* temp = xml.element("temporal", 0);
695  m_temporal = xml_temporal(*temp);
696  }
697  else {
699  m_temporal = temporal.clone();
700  }
701 
702  // Read model attributes
703  read_attributes(xml);
704 
705  // Set parameter pointers
706  set_pointers();
707 
708  // Return
709  return;
710 }
711 
712 
713 /***********************************************************************//**
714  * @brief Write model into XML element
715  *
716  * @param[in] xml XML element.
717  *
718  * Writes the sky model into an XML source library. The format written to
719  * the @p xml element is as follows:
720  *
721  * <source name=".." type=".." instrument=".." id="..">
722  * <spectrum type="..">
723  * ..
724  * </spectrum>
725  * <spatialModel type="..">
726  * ..
727  * </spatialModel>
728  * <temporal type="..">
729  * ..
730  * </temporal>
731  * </source>
732  *
733  * For compatibility reasons the temporal element will only be written if it
734  * is a non-constant component or a constant component with a normalization
735  * that differs from unity.
736  ***************************************************************************/
738 {
739  // Initialise pointer on source
740  GXmlElement* src = NULL;
741 
742  // Search corresponding source
743  int n = xml.elements("source");
744  for (int k = 0; k < n; ++k) {
745  GXmlElement* element = xml.element("source", k);
746  if (element->attribute("name") == name()) {
747  src = element;
748  break;
749  }
750  }
751 
752  // If the temporal model is not a constant with unit normalization then
753  // set cons to a NULL pointer
754  GModelTemporalConst* cons = dynamic_cast<GModelTemporalConst*>(temporal());
755  if (cons != NULL) {
756  if (cons->norm() != 1.0) {
757  cons = NULL;
758  }
759  }
760 
761  // If no source with corresponding name was found then append one
762  if (src == NULL) {
763  src = xml.append("source");
764  if (spatial() != NULL) src->append(GXmlElement("spatialModel"));
765  if (spectral() != NULL) src->append(GXmlElement("spectrum"));
766  if (temporal() != NULL && cons == NULL) src->append(GXmlElement("temporal"));
767  }
768 
769  // Write spatial model
770  if (spatial() != NULL) {
771  GXmlElement* spat = src->element("spatialModel", 0);
772  spatial()->write(*spat);
773  }
774 
775  // Write spectral model
776  if (spectral() != NULL) {
777  GXmlElement* spec = src->element("spectrum", 0);
778  spectral()->write(*spec);
779  }
780 
781  // Write temporal model (only if not a constant with unit normalization
782  // factor)
783  if (temporal() != NULL && cons == NULL) {
784  GXmlElement* temp = src->element("temporal", 0);
785  temporal()->write(*temp);
786  }
787 
788  // Write model attributes
789  write_attributes(*src);
790 
791  // Return
792  return;
793 }
794 
795 
796 /***********************************************************************//**
797  * @brief Print model information
798  *
799  * @param[in] chatter Chattiness.
800  * @return String containing model information.
801  ***************************************************************************/
802 std::string GCTAModelBackground::print(const GChatter& chatter) const
803 {
804  // Initialise result string
805  std::string result;
806 
807  // Continue only if chatter is not silent
808  if (chatter != SILENT) {
809 
810  // Append header
811  result.append("=== GCTAModelBackground ===");
812 
813  // Determine number of parameters per type
814  int n_spatial = (spatial() != NULL) ? spatial()->size() : 0;
815  int n_spectral = (spectral() != NULL) ? spectral()->size() : 0;
816  int n_temporal = (temporal() != NULL) ? temporal()->size() : 0;
817 
818  // Append attributes
819  result.append("\n"+print_attributes());
820 
821  // Append model type
822  result.append("\n"+gammalib::parformat("Model type")+type());
823 
824  // Append model type
825  result.append("\n"+gammalib::parformat("Model components"));
826  if (n_spatial > 0) {
827  result.append("\""+spatial()->type()+"\"");
828  if (n_spectral > 0 || n_temporal > 0) {
829  result.append(" * ");
830  }
831  }
832  if (n_spectral > 0) {
833  result.append("\""+spectral()->type()+"\"");
834  if (n_temporal > 0) {
835  result.append(" * ");
836  }
837  }
838  if (n_temporal > 0) {
839  result.append("\""+temporal()->type()+"\"");
840  }
841 
842  // Append parameters
843  result.append("\n"+gammalib::parformat("Number of parameters") +
844  gammalib::str(size()));
845  result.append("\n"+gammalib::parformat("Number of spatial par's") +
846  gammalib::str(n_spatial));
847  for (int i = 0; i < n_spatial; ++i) {
848  result.append("\n"+(*spatial())[i].print());
849  }
850  result.append("\n"+gammalib::parformat("Number of spectral par's") +
851  gammalib::str(n_spectral));
852  for (int i = 0; i < n_spectral; ++i) {
853  result.append("\n"+(*spectral())[i].print());
854  }
855  result.append("\n"+gammalib::parformat("Number of temporal par's") +
856  gammalib::str(n_temporal));
857  for (int i = 0; i < n_temporal; ++i) {
858  result.append("\n"+(*temporal())[i].print());
859  }
860 
861  } // endif: chatter was not silent
862 
863  // Return result
864  return result;
865 }
866 
867 
868 /*==========================================================================
869  = =
870  = Private methods =
871  = =
872  ==========================================================================*/
873 
874 /***********************************************************************//**
875  * @brief Initialise class members
876  ***************************************************************************/
878 {
879  // Initialise members
880  m_spatial = NULL;
881  m_spectral = NULL;
882  m_temporal = NULL;
883 
884  // Return
885  return;
886 }
887 
888 
889 /***********************************************************************//**
890  * @brief Copy class members
891  *
892  * @param[in] model Background model.
893  ***************************************************************************/
895 {
896  // Clone radial, spectral and temporal model components
897  m_spatial = (model.m_spatial != NULL) ? model.m_spatial->clone() : NULL;
898  m_spectral = (model.m_spectral != NULL) ? model.m_spectral->clone() : NULL;
899  m_temporal = (model.m_temporal != NULL) ? model.m_temporal->clone() : NULL;
900 
901  // Set parameter pointers
902  set_pointers();
903 
904  // Return
905  return;
906 }
907 
908 
909 /***********************************************************************//**
910  * @brief Delete class members
911  ***************************************************************************/
913 {
914  // Free memory
915  if (m_spatial != NULL) delete m_spatial;
916  if (m_spectral != NULL) delete m_spectral;
917  if (m_temporal != NULL) delete m_temporal;
918 
919  // Signal free pointers
920  m_spatial = NULL;
921  m_spectral = NULL;
922  m_temporal = NULL;
923 
924  // Return
925  return;
926 }
927 
928 
929 /***********************************************************************//**
930  * @brief Set pointers
931  *
932  * Set pointers to all model parameters. The pointers are stored in a vector
933  * that is member of the GModelData base class.
934  ***************************************************************************/
936 {
937  // Clear parameters
938  m_pars.clear();
939 
940  // Determine the number of parameters
941  int n_spatial = (spatial() != NULL) ? spatial()->size() : 0;
942  int n_spectral = (spectral() != NULL) ? spectral()->size() : 0;
943  int n_temporal = (temporal() != NULL) ? temporal()->size() : 0;
944  int n_pars = n_spatial + n_spectral + n_temporal;
945 
946  // Continue only if there are parameters
947  if (n_pars > 0) {
948 
949  // Gather spatial parameter pointers
950  for (int i = 0; i < n_spatial; ++i) {
951  m_pars.push_back(&((*spatial())[i]));
952  }
953 
954  // Gather spectral parameters
955  for (int i = 0; i < n_spectral; ++i) {
956  m_pars.push_back(&((*spectral())[i]));
957  }
958 
959  // Gather temporal parameters
960  for (int i = 0; i < n_temporal; ++i) {
961  m_pars.push_back(&((*temporal())[i]));
962  }
963 
964  }
965 
966  // Return
967  return;
968 }
969 
970 
971 /***********************************************************************//**
972  * @brief Return pointer to spatial model from XML element
973  *
974  * @param[in] spatial XML element.
975  * @return Pointer to spatial model.
976  *
977  * Returns pointer to a spatial model that is defined in an XML element.
978  ***************************************************************************/
980 {
981  // Get radial model
982  GCTAModelSpatialRegistry registry;
983  GCTAModelSpatial* ptr = registry.alloc(spatial);
984 
985  // Return pointer
986  return ptr;
987 }
988 
989 
990 /***********************************************************************//**
991  * @brief Return pointer to spectral model from XML element
992  *
993  * @param[in] spectral XML element.
994  * @return Pointer to spectral model.
995  *
996  * Returns pointer to spectral model that is defined in an XML element.
997  ***************************************************************************/
999 {
1000  // Get spectral model
1001  GModelSpectralRegistry registry;
1002  GModelSpectral* ptr = registry.alloc(spectral);
1003 
1004  // Return pointer
1005  return ptr;
1006 }
1007 
1008 
1009 /***********************************************************************//**
1010  * @brief Return pointer to temporal model from XML element
1011  *
1012  * @param[in] temporal XML element.
1013  * @return Pointer to temporal model.
1014  *
1015  * Returns pointer to temporal model that is defined in an XML element.
1016  ***************************************************************************/
1018 {
1019  // Get temporal model
1020  GModelTemporalRegistry registry;
1021  GModelTemporal* ptr = registry.alloc(temporal);
1022 
1023  // Return pointer
1024  return ptr;
1025 }
1026 
1027 
1028 /***********************************************************************//**
1029  * @brief Verifies if model has all components
1030  *
1031  * Returns 'true' if models has a spatial, a spectral and a temporal
1032  * component. Otherwise returns 'false'.
1033  ***************************************************************************/
1035 {
1036  // Set result
1037  bool result = ((spatial() != NULL) &&
1038  (spectral() != NULL) &&
1039  (temporal() != NULL));
1040 
1041  // Return result
1042  return result;
1043 }
virtual void write(GXmlElement &xml) const =0
Constant temporal model class interface definition.
GCTAModelBackground(void)
Void constructor.
virtual GCTAInstDir mc(const GEnergy &energy, const GTime &time, const GCTAObservation &obs, GRan &ran) const
Returns MC instrument direction.
#define G_MC
Abstract spatial model class.
void free_members(void)
Delete class members.
Definition: GModel.cpp:660
GModelTemporal * xml_temporal(const GXmlElement &temporal) const
Return pointer to temporal model from XML element.
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
const GCTAInstDir & cta_dir(const std::string &origin, const GEvent &event)
Retrieve CTA instrument direction from generic event.
int size(void) const
Return number of times.
Definition: GTimes.hpp:100
std::string print_attributes(void) const
Print model attributes.
Definition: GModel.cpp:752
Abstract spectral model base class.
Interface definition for the model registry class.
virtual void read(const GXmlElement &xml)
Read model from XML element.
Abstract temporal model base class.
virtual void roi(const GRoi &roi)
Set ROI.
GModelTemporal * m_temporal
Temporal model.
int size(void) const
Return number of parameters.
void ebounds(const GEbounds &ebounds)
Set energy boundaries.
Definition: GEvents.cpp:136
int size(void) const
Return number of energy boundaries.
Definition: GEbounds.hpp:157
virtual double eval(const GEnergy &srcEng, const GTime &srcTime=GTime(), const bool &gradients=false) const =0
virtual void write(GXmlElement &xml) const =0
void append(const GCTAEventAtom &event)
Append event to event list.
virtual double npred(const GEnergy &energy, const GTime &time, const GObservation &obs) const
Return spatially integrated background rate in units of events MeV s .
const GCTAPointing & cta_pnt(const std::string &origin, const GObservation &obs)
Retrieve CTA pointing from generic observation.
bool valid_model(void) const
Verifies if model has all components.
Abstract interface for the event classes.
Definition: GEvent.hpp:71
CTA event list class.
XML element node class.
Definition: GXmlElement.hpp:47
void gti(const GGti &gti)
Set Good Time Intervals.
Definition: GEvents.cpp:154
CTA pointing class interface definition.
void free_members(void)
Delete class members.
Definition: GModelData.cpp:181
Spectral model registry class definition.
virtual void write(GXmlElement &xml) const =0
GModelSpectral * spectral(void) const
Return spectral model component.
Random number generator class.
Definition: GRan.hpp:44
Interface for the CTA region of interest class.
Definition: GCTARoi.hpp:48
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.
virtual const GTime & time(void) const =0
GModelSpectral * m_spectral
Spectral model.
int size(void) const
Return number of Good Time Intervals.
Definition: GGti.hpp:153
Interface definition for the spatial model registry class.
Time container class.
Definition: GTimes.hpp:45
std::vector< GModelPar * > m_pars
Pointers to all model parameters.
Definition: GModel.hpp:171
int size(void) const
Return number of parameters in model.
Definition: GModel.hpp:225
void init_members(void)
Initialise class members.
Definition: GModel.cpp:612
void copy_members(const GCTAModelBackground &model)
Copy class members.
Definition of support function used by CTA classes.
void free_members(void)
Delete class members.
virtual GCTAEventList * mc(const GObservation &obs, GRan &ran) const
Return simulated list of events.
Energy boundaries container class.
Definition: GEbounds.hpp:60
const GXmlAttribute * attribute(const int &index) const
Return attribute.
const std::string & name(void) const
Return parameter name.
Definition: GModel.hpp:253
void reserve(const int &number)
Reserves space for events.
const GEnergy & emin(void) const
Return minimum energy of all intervals.
Definition: GEbounds.hpp:181
GCTAModelSpatial * spatial(void) const
Return spatial model component.
void set_pointers(void)
Set pointers.
int size(void) const
Return number of parameters.
GModelTemporal * alloc(const GXmlElement &xml) const
Allocate temporal model that is found in XML element.
virtual const GEnergy & energy(void) const =0
CTA pointing class.
Abstract data model class.
Definition: GModelData.hpp:53
CTA instrument direction class interface definition.
const GCTAEventList & cta_event_list(const std::string &origin, const GObservation &obs)
Retrieve CTA event list from generic observation.
GModelSpectral * alloc(const GXmlElement &xml) const
Allocate spectral model that is found in XML element.
virtual GCTAModelBackground & operator=(const GCTAModelBackground &model)
Assignment operator.
virtual GModelData & operator=(const GModelData &model)
Assignment operator.
Definition: GModelData.cpp:120
int size(void) const
Return number of model parameters.
GChatter
Definition: GTypemaps.hpp:33
Good Time Interval class.
Definition: GGti.hpp:62
GCTAModelSpatial * xml_spatial(const GXmlElement &spatial) const
Return pointer to spatial model from XML element.
virtual double eval(const GTime &srcTime, const bool &gradients=false) const =0
const GCTAObservation & cta_obs(const std::string &origin, const GObservation &obs)
Retrieve CTA observation from generic observation.
Abstract observation base class.
void init_members(void)
Initialise class members.
const GTime & tstop(void) const
Returns latest stop time in Good Time Intervals.
Definition: GGti.hpp:206
const GTime & tstart(void) const
Returns earliest start time in Good Time Intervals.
Definition: GGti.hpp:193
CTA region of interest class interface definition.
CTA observation class interface definition.
CTA exception handler interface definition.
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:708
virtual GModelSpectral * clone(void) const =0
Clones object.
Background model class.
virtual GCTAModelBackground * clone(void) const
Clone instance.
virtual void write(GXmlElement &xml) const
Write model into XML element.
virtual GXmlElement * element(const int &index)
Return pointer to GXMLElement child.
Definition: GXmlNode.cpp:634
GCTAModelSpatial * m_spatial
Spatial model.
GModelTemporal * temporal(void) const
Return temporal model component.
Temporal model registry class definition.
const GCTAModelBackground g_cta_background_seed
CTA event atom class.
Exception handler interface definition.
virtual double eval(const GCTAInstDir &dir, const GEnergy &energy, const GTime &time, const bool &gradients=false) const =0
virtual GModelTemporal * clone(void) const =0
Clones object.
void read_attributes(const GXmlElement &xml)
Read model attributes.
Definition: GModel.cpp:672
CTA instrument direction class.
Definition: GCTAInstDir.hpp:63
GCTAModelSpatial * alloc(const GXmlElement &xml) const
Allocate spatial model that is found in XML element.
double norm(void) const
Return normalization factor.
const GEnergy & emax(void) const
Return maximum energy of all intervals.
Definition: GEbounds.hpp:193
virtual double flux(const GEnergy &emin, const GEnergy &emax) const =0
Model registry class definition.
virtual GCTAModelSpatial * clone(void) const =0
Clones object.
virtual GXmlNode * append(const GXmlNode &node)
Append XML child node.
Definition: GXmlNode.cpp:284
CTA observation class.
std::string parformat(const std::string &s, const int &indent=0)
Convert string in parameter format.
Definition: GTools.cpp:1033
#define G_EVAL
GModelSpectral * xml_spectral(const GXmlElement &spectral) const
Return pointer to spectral model from XML element.
virtual double npred(const GEnergy &energy, const GTime &time, const GObservation &obs) const
Return integral of spatial model component.
Integration class interface definition.
virtual double eval(const GEvent &event, const GObservation &obs, const bool &gradients=false) const
Return background rate in units of events MeV s sr .
Spatial model registry class definition.
virtual std::string print(const GChatter &chatter=NORMAL) const
Print model information.
void init_members(void)
Initialise class members.
Definition: GModelData.cpp:159
virtual ~GCTAModelBackground(void)
Destructor.
Constant temporal model class.
const GCTAInstDir & dir(void) const
Return instrument direction.
Mathematical function definitions.
virtual void clear(void)
Clear instance.
Class that handles energies in a unit independent way.
Definition: GEnergy.hpp:48
virtual std::string type(void) const
Return model type.
Background model class interface definition.
std::string str(const unsigned short int &value)
Convert unsigned short integer value into string.
Definition: GTools.cpp:415