GammaLib  2.1.0.dev
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
GCTAModelAeffBackground.cpp
Go to the documentation of this file.
1 /***************************************************************************
2  * GCTAModelAeffBackground.cpp - CTA Aeff background model class *
3  * ----------------------------------------------------------------------- *
4  * copyright (C) 2015-2020 by Michael Mayer *
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 GCTAModelAeffBackground.cpp
23  * @brief CTA Aeff background model class implementation
24  * @author Michael Mayer
25  */
26 
27 /* __ Includes ___________________________________________________________ */
28 #ifdef HAVE_CONFIG_H
29 #include <config.h>
30 #endif
31 #include "GTools.hpp"
32 #include "GMath.hpp"
33 #include "GIntegral.hpp"
34 #include "GModelRegistry.hpp"
37 #include "GModelTemporalConst.hpp"
38 #include "GModelSpectralNodes.hpp"
40 #include "GCTAObservation.hpp"
41 #include "GCTAResponseIrf.hpp"
42 #include "GCTABackground.hpp"
43 #include "GCTASupport.hpp"
44 
45 /* __ Globals ____________________________________________________________ */
47 const GModelRegistry g_cta_aeff_background_registry(&g_cta_aeff_background_seed);
48 
49 /* __ Method name definitions ____________________________________________ */
50 #define G_EVAL "GCTAModelAeffBackground::eval(GEvent&, GObservation&, bool&)"
51 #define G_NPRED "GCTAModelAeffBackground::npred(GEnergy&, GTime&,"\
52  " GObservation&)"
53 #define G_MC "GCTAModelAeffBackground::mc(GObservation&, GRan&)"
54 #define G_XML_SPECTRAL "GCTAModelAeffBackground::xml_spectral(GXmlElement&)"
55 #define G_XML_TEMPORAL "GCTAModelAeffBackground::xml_temporal(GXmlElement&)"
56 #define G_AEFF_INTEGRAL "GCTAModelAeffBackground::aeff_integral("\
57  "GObservation&, double&)"
58 
59 /* __ Macros _____________________________________________________________ */
60 
61 /* __ Coding definitions _________________________________________________ */
62 #define G_USE_NPRED_CACHE
63 
64 /* __ Debug definitions __________________________________________________ */
65 //#define G_DUMP_MC
66 //#define G_DEBUG_NPRED
67 
68 /* __ Constants __________________________________________________________ */
69 
70 
71 /*==========================================================================
72  = =
73  = Constructors/destructors =
74  = =
75  ==========================================================================*/
76 
77 /***********************************************************************//**
78  * @brief Void constructor
79  ***************************************************************************/
81 {
82  // Initialise class members
83  init_members();
84 
85  // Return
86  return;
87 }
88 
89 
90 /***********************************************************************//**
91  * @brief XML constructor
92  *
93  * @param[in] xml XML element.
94  *
95  * Constructs a CTA effective area background model from the information
96  * provided by an XML element. The XML element is expected to have the
97  * following structure
98  *
99  * <source name="..." type="..." instrument="...">
100  * <spectrum type="...">
101  * ...
102  * </spectrum>
103  * </source>
104  *
105  * Optionally, a temporal model may be provided using the following
106  * syntax
107  *
108  * <source name="..." type="..." instrument="...">
109  * <spectrum type="...">
110  * ...
111  * </spectrum>
112  * <temporalModel type="...">
113  * ...
114  * </temporalModel>
115  * </source>
116  *
117  * If no temporal component is found a constant model is assumed.
118  ***************************************************************************/
120  GModelData(xml)
121 {
122  // Initialise members
123  init_members();
124 
125  // Read XML
126  read(xml);
127 
128  // Set parameter pointers
129  set_pointers();
130 
131  // Return
132  return;
133 }
134 
135 
136 /***********************************************************************//**
137  * @brief Construct from spectral component
138  *
139  * @param[in] spectral Spectral model component.
140  *
141  * Constructs a CTA effective area background model from a @p spectral
142  * model component. The temporal component is assumed to be constant.
143  ***************************************************************************/
145  GModelData()
146 {
147  // Initialise members
148  init_members();
149 
150  // Allocate temporal constant model
152 
153  // Clone model components
154  m_spectral = spectral.clone();
155  m_temporal = temporal.clone();
156 
157  // Set parameter pointers
158  set_pointers();
159 
160  // Return
161  return;
162 }
163 
164 
165 /***********************************************************************//**
166  * @brief Construct from model components
167  *
168  * @param[in] spectral Spectral model component.
169  * @param[in] temporal Temporal model component.
170  *
171  * Constructs a CTA effective area background model from a @p spectral and a
172  * @p temporal component.
173  ***************************************************************************/
175  const GModelTemporal& temporal) :
176  GModelData()
177 {
178  // Initialise members
179  init_members();
180 
181  // Clone model components
182  m_spectral = spectral.clone();
183  m_temporal = temporal.clone();
184 
185  // Set parameter pointers
186  set_pointers();
187 
188  // Return
189  return;
190 }
191 
192 
193 /***********************************************************************//**
194  * @brief Copy constructor
195  *
196  * @param[in] aeff CTA instrument background model.
197  ***************************************************************************/
199  GModelData(aeff)
200 {
201  // Initialise class members
202  init_members();
203 
204  // Copy members
205  copy_members(aeff);
206 
207  // Return
208  return;
209 }
210 
211 
212 /***********************************************************************//**
213  * @brief Destructor
214  ***************************************************************************/
216 {
217  // Free members
218  free_members();
219 
220  // Return
221  return;
222 }
223 
224 
225 /*==========================================================================
226  = =
227  = Operators =
228  = =
229  ==========================================================================*/
230 
231 /***********************************************************************//**
232  * @brief Assignment operator
233  *
234  * @param[in] aeff CTA effective area background model.
235  * @return CTA effective area background model.
236  ***************************************************************************/
238 {
239  // Execute only if object is not identical
240  if (this != &aeff) {
241 
242  // Copy base class members
243  this->GModelData::operator=(aeff);
244 
245  // Free members
246  free_members();
247 
248  // Initialise private members
249  init_members();
250 
251  // Copy members
252  copy_members(aeff);
253 
254  } // endif: object was not identical
255 
256  // Return this object
257  return *this;
258 }
259 
260 
261 /*==========================================================================
262  = =
263  = Public methods =
264  = =
265  ==========================================================================*/
266 
267 /***********************************************************************//**
268  * @brief Clear CTA effective area background model
269  *
270  * This method properly resets the CTA effective area background model to an
271  * initial state.
272  ***************************************************************************/
274 {
275  // Free class members (base and derived classes, derived class first)
276  free_members();
277  this->GModelData::free_members();
278 
279  // Initialise members
280  this->GModelData::init_members();
281  init_members();
282 
283  // Return
284  return;
285 }
286 
287 
288 /***********************************************************************//**
289  * @brief Clone CTA effective area background model
290  *
291  * @return Pointer to deep copy of CTA effective area background model.
292  ***************************************************************************/
294 {
295  return new GCTAModelAeffBackground(*this);
296 }
297 
298 
299 /***********************************************************************//**
300  * @brief Set spectral model component
301  *
302  * @param[in] spectral Pointer to spectral model component.
303  *
304  * Sets the spectral model component of the model.
305  ***************************************************************************/
307 {
308  // Free spectral model component only if it differs from current
309  // component. This prevents unintential deallocation of the argument
310  if ((m_spectral != NULL) && (m_spectral != spectral)) {
311  delete m_spectral;
312  }
313 
314  // Clone spectral model component if it exists, otherwise set pointer
315  // to NULL
316  m_spectral = (spectral != NULL) ? spectral->clone() : NULL;
317 
318  // Set parameter pointers
319  set_pointers();
320 
321  // Return
322  return;
323 }
324 
325 
326 /***********************************************************************//**
327  * @brief Set temporal model component
328  *
329  * @param[in] temporal Pointer to temporal model component.
330  *
331  * Sets the temporal model component of the model.
332  ***************************************************************************/
334 {
335  // Free temporal model component only if it differs from current
336  // component. This prevents unintential deallocation of the argument
337  if ((m_temporal != NULL) && (m_temporal != temporal)) {
338  delete m_temporal;
339  }
340 
341  // Clone temporal model component if it exists, otherwise set pointer
342  // to NULL
343  m_temporal = (temporal != NULL) ? temporal->clone() : NULL;
344 
345  // Set parameter pointers
346  set_pointers();
347 
348  // Return
349  return;
350 }
351 
352 
353 /***********************************************************************//**
354  * @brief Return background rate in units of events MeV\f$^{-1}\f$
355  * s\f$^{-1}\f$ sr\f$^{-1}\f$
356  *
357  * @param[in] event Observed event.
358  * @param[in] obs Observation.
359  * @param[in] gradients Compute gradients?
360  * @return Background rate (events MeV\f$^{-1}\f$ s\f$^{-1}\f$ sr\f$^{-1}\f$).
361  *
362  * The method returns a real rate, defined by the number of counts per MeV,
363  * steradian and ontime.
364  *
365  * If the @p gradients flag is true the method will also set the parameter
366  * gradients of the model parameters.
367  ***************************************************************************/
369  const GObservation& obs,
370  const bool& gradients) const
371 {
372  // Get reference on CTA pointing and effective area response from
373  // observation and reference on CTA instrument direction from event
374  const GCTAPointing& pnt = gammalib::cta_pnt(G_EVAL, obs);
375  const GCTAAeff& aeff = gammalib::cta_rsp_aeff(G_EVAL, obs);
376  const GCTAInstDir& dir = gammalib::cta_dir(G_EVAL, event);
377 
378  // Set instrument direction
379  GCTAInstDir inst_dir = pnt.instdir(dir.dir());
380 
381  // Evaluate function
382  double logE = event.energy().log10TeV();
383  double spat = aeff(logE,
384  inst_dir.theta(),
385  inst_dir.phi(),
386  pnt.zenith(),
387  pnt.azimuth(), false);
388  double spec = (spectral() != NULL)
389  ? spectral()->eval(event.energy(), event.time(), gradients)
390  : 1.0;
391  double temp = (temporal() != NULL)
392  ? temporal()->eval(event.time(), gradients) : 1.0;
393 
394  // Compute value
395  double value = spat * spec * temp;
396 
397  // Optionally compute partial derivatives
398  if (gradients) {
399 
400  // Multiply factors to spectral gradients
401  if (spectral() != NULL) {
402  double fact = spat * temp;
403  if (fact != 1.0) {
404  for (int i = 0; i < spectral()->size(); ++i)
405  (*spectral())[i].factor_gradient((*spectral())[i].factor_gradient() * fact );
406  }
407  }
408 
409  // Multiply factors to temporal gradients
410  if (temporal() != NULL) {
411  double fact = spat * spec;
412  if (fact != 1.0) {
413  for (int i = 0; i < temporal()->size(); ++i)
414  (*temporal())[i].factor_gradient((*temporal())[i].factor_gradient() * fact );
415  }
416  }
417 
418  } // endif: computed partial derivatives
419 
420  // Return value
421  return value;
422 }
423 
424 
425 /***********************************************************************//**
426  * @brief Return spatially integrated background rate in units of
427  * events MeV\f$^{-1}\f$ s\f$^{-1}\f$
428  *
429  * @param[in] obsEng Measured event energy.
430  * @param[in] obsTime Measured event time.
431  * @param[in] obs Observation.
432  * @return Spatially integrated background rate
433  * (events MeV\f$^{-1}\f$ s\f$^{-1}\f$)
434  *
435  * Spatially integrates the effective area background model for a given
436  * measured event energy and event time. The method returns a real rate,
437  * defined as the number of counts per MeV and ontime.
438  ***************************************************************************/
440  const GTime& obsTime,
441  const GObservation& obs) const
442 {
443  // Initialise result
444  double npred = 0.0;
445  bool has_npred = false;
446 
447  // Build unique identifier
448  std::string id = obs.instrument() + "::" + obs.id();
449 
450  // Check if Npred value is already in cache
451  #if defined(G_USE_NPRED_CACHE)
452  if (!m_npred_names.empty()) {
453 
454  // Search for unique identifier, and if found, recover Npred value
455  // and break
456  for (int i = 0; i < m_npred_names.size(); ++i) {
457  if (m_npred_names[i] == id && m_npred_energies[i] == obsEng) {
458  npred = m_npred_values[i];
459  has_npred = true;
460  #if defined(G_DEBUG_NPRED)
461  std::cout << "GCTAModelAeffBackground::npred:";
462  std::cout << " cache=" << i;
463  std::cout << " npred=" << npred << std::endl;
464  #endif
465  break;
466  }
467  }
468 
469  } // endif: there were values in the Npred cache
470  #endif
471 
472  // Continue only if no Npred cache value has been found
473  if (!has_npred) {
474 
475  // Evaluate only if model is valid
476  if (valid_model()) {
477 
478  // Get log10 of energy in TeV
479  double logE = obsEng.log10TeV();
480 
481  // Spatially integrate effective area component
482  npred = this->aeff_integral(obs, logE);
483 
484  // Store result in Npred cache
485  #if defined(G_USE_NPRED_CACHE)
486  m_npred_names.push_back(id);
487  m_npred_energies.push_back(obsEng);
488  m_npred_times.push_back(obsTime);
489  m_npred_values.push_back(npred);
490  #endif
491 
492  // Debug: Check for NaN
493  #if defined(G_NAN_CHECK)
494  if (gammalib::is_notanumber(npred) || gammalib::is_infinite(npred)) {
495  std::string origin = "GCTAModelAeffBackground::npred";
496  std::string message = " NaN/Inf encountered (npred=" +
497  gammalib::str(npred) + ")";
498  gammalib::warning(origin, message);
499  }
500  #endif
501 
502  } // endif: model was valid
503 
504  } // endif: Npred computation required
505 
506  // Multiply in spectral and temporal components
507  npred *= spectral()->eval(obsEng, obsTime);
508  npred *= temporal()->eval(obsTime);
509 
510  // Return Npred
511  return npred;
512 }
513 
514 
515 /***********************************************************************//**
516  * @brief Return simulated list of events
517  *
518  * @param[in] obs Observation.
519  * @param[in] ran Random number generator.
520  * @return Pointer to list of simulated events (needs to be de-allocated by
521  * client)
522  *
523  * Draws a sample of events from the background model using a Monte
524  * Carlo simulation. The region of interest, the energy boundaries and the
525  * good time interval for the sampling will be extracted from the observation
526  * argument that is passed to the method. The method also requires a random
527  * number generator of type GRan which is passed by reference, hence the
528  * state of the random number generator will be changed by the method.
529  *
530  * For each event in the returned event list, the sky direction, the nominal
531  * coordinates (DETX and DETY), the energy and the time will be set.
532  ***************************************************************************/
534  GRan& ran) const
535 {
536  // Initialise new event list
537  GCTAEventList* list = new GCTAEventList;
538 
539  // Continue only if model is valid)
540  if (valid_model()) {
541 
542  // Retrieve CTA pointing, effective area response and event list
543  const GCTAPointing& pnt = gammalib::cta_pnt(G_MC, obs);
544  const GCTAAeff& aeff = gammalib::cta_rsp_aeff(G_MC, obs);
545  const GCTAEventList& events = gammalib::cta_event_list(G_MC, obs);
546 
547  // Get simulation region
548  const GCTARoi& roi = events.roi();
549  const GEbounds& ebounds = events.ebounds();
550  const GGti& gti = events.gti();
551 
552  // Get maximum offset value for simulations
553  double max_theta = pnt.dir().dist(roi.centre().dir()) +
554  roi.radius() * gammalib::deg2rad;
555  double cos_max_theta = std::cos(max_theta);
556 
557  // Set simulation region for result event list
558  list->roi(roi);
559  list->ebounds(ebounds);
560  list->gti(gti);
561 
562  // Set up spectral model to draw random energies from. Here we use
563  // a fixed energy sampling for an instance of GModelSpectralNodes.
564  // This is analogous to to the GCTAModelIrfBackground::mc method.
565  // We make sure that only non-negative nodes get appended.
566  GEbounds spectral_ebounds(m_n_mc_energies,
567  ebounds.emin(),
568  ebounds.emax(),
569  "LOG");
571  for (int i = 0; i < spectral_ebounds.size(); ++i) {
572  GEnergy energy = spectral_ebounds.elogmean(i);
573  double intensity = aeff_integral(obs, energy.log10TeV());
574  double norm = m_spectral->eval(energy, events.tstart());
575  double arg = norm * intensity;
576  if (arg > 0.0) {
577  spectral.append(energy, arg);
578  }
579  }
580 
581  // Loop over all energy bins
582  for (int ieng = 0; ieng < ebounds.size(); ++ieng) {
583 
584  // Compute the background rate in model within the energy
585  // boundaries from spectral component (units: cts/s).
586  // Note that the time here is ontime. Deadtime correction will
587  // be done later.
588  double rate = spectral.flux(ebounds.emin(ieng), ebounds.emax(ieng));
589 
590  // Debug option: dump rate
591  #if defined(G_DUMP_MC)
592  std::cout << "GCTAModelAeffBackground::mc(\"" << name() << "\": ";
593  std::cout << "rate=" << rate << " cts/s)" << std::endl;
594  #endif
595 
596  // If the rate is not positive then skip this energy bins
597  if (rate <= 0.0) {
598  continue;
599  }
600 
601  // Loop over all good time intervals
602  for (int itime = 0; itime < gti.size(); ++itime) {
603 
604  // Get Monte Carlo event arrival times from temporal model
605  GTimes times = m_temporal->mc(rate,
606  gti.tstart(itime),
607  gti.tstop(itime),
608  ran);
609 
610  // Get number of events
611  int n_events = times.size();
612 
613  // Reserve space for events
614  if (n_events > 0) {
615  list->reserve(n_events);
616  }
617 
618  // Debug option: provide number of times and initialize
619  // statisics
620  #if defined(G_DUMP_MC)
621  std::cout << " Interval " << itime;
622  std::cout << " times=" << n_events << std::endl;
623  int n_killed_by_roi = 0;
624  #endif
625 
626  // Loop over events
627  for (int i = 0; i < n_events; ++i) {
628 
629  // Get Monte Carlo event energy from spectral model
630  GEnergy energy = spectral.mc(ebounds.emin(ieng),
631  ebounds.emax(ieng),
632  times[i],
633  ran);
634 
635  // Get maximum effective area for rejection method
636  double max_aeff = aeff.max(energy.log10TeV(), pnt.zenith(),
637  pnt.azimuth(), false);
638 
639  // Skip event if the maximum effective area is not positive
640  if (max_aeff <= 0.0) {
641  continue;
642  }
643 
644  // Initialise randomised coordinates
645  double offset = 0.0;
646  double phi = 0.0;
647 
648  // Initialise acceptance fraction and counter of zeros for
649  // rejection method
650  double acceptance_fraction = 0.0;
651  int zeros = 0;
652 
653  // Start rejection method loop
654  do {
655 
656  // Throw random offset and azimuth angle in
657  // considered range
658  offset = std::acos(1.0 - ran.uniform() *
659  (1.0 - cos_max_theta));
660  phi = ran.uniform() * gammalib::twopi;
661 
662  // Compute function value at this offset angle
663  double value = aeff(energy.log10TeV(), offset, phi,
664  pnt.zenith(), pnt.azimuth(),
665  false);
666 
667  // If the value is not positive then increment the
668  // zeros counter and fall through. The counter assures
669  // that this loop does not lock up.
670  if (value <= 0.0) {
671  zeros++;
672  continue;
673  }
674 
675  // Value is non-zero so reset the zeros counter
676  zeros = 0;
677 
678  // Compute acceptance fraction
679  acceptance_fraction = value / max_aeff;
680 
681  } while ((ran.uniform() > acceptance_fraction) &&
682  (zeros < 1000));
683 
684  // If the zeros counter is non-zero then the loop was
685  // exited due to exhaustion and the event is skipped
686  if (zeros > 0) {
687  continue;
688  }
689 
690  // Rotate pointing direction by offset and azimuth angle
691  GSkyDir skydir = pnt.dir();
692  skydir.rotate_deg(phi * gammalib::rad2deg,
693  offset * gammalib::rad2deg);
694 
695  // Convert rotated pointing direction in instrument system
696  GCTAInstDir mc_dir = pnt.instdir(skydir);
697 
698  // Allocate event
699  GCTAEventAtom event;
700 
701  // Set event attributes
702  event.dir(mc_dir);
703  event.energy(energy);
704  event.time(times[i]);
705 
706  // Append event to list if it falls in ROI
707  if (events.roi().contains(event)) {
708  list->append(event);
709  }
710  #if defined(G_DUMP_MC)
711  else {
712  n_killed_by_roi++;
713  }
714  #endif
715 
716  } // endfor: looped over all events
717 
718  // Debug option: provide statisics
719  #if defined(G_DUMP_MC)
720  std::cout << " Killed by ROI=";
721  std::cout << n_killed_by_roi << std::endl;
722  #endif
723 
724  } // endfor: looped over all GTIs
725 
726  } // endfor: looped over all energy boundaries
727 
728  } // endif: model was valid
729 
730  // Return
731  return list;
732 }
733 
734 
735 /***********************************************************************//**
736  * @brief Read CTA effective area background model from XML element
737  *
738  * @param[in] xml XML element.
739  *
740  * Set up CTA effective area background model from the information provided by
741  * an XML element. The XML element is expected to have the following
742  * structure
743  *
744  * <source name="..." type="..." instrument="...">
745  * <spectrum type="...">
746  * ...
747  * </spectrum>
748  * </source>
749  *
750  * Optionally, a temporal model may be provided using the following
751  * syntax
752  *
753  * <source name="..." type="..." instrument="...">
754  * <spectrum type="...">
755  * ...
756  * </spectrum>
757  * <temporal type="...">
758  * ...
759  * </temporal>
760  * </source>
761  *
762  * If no temporal component is found a constant model is assumed.
763  ***************************************************************************/
765 {
766  // Clear model
767  clear();
768 
769  // Initialise XML elements
770  const GXmlElement* spectral = NULL;
771 
772  // Get pointer on spectrum
773  spectral = xml.element("spectrum", 0);
774 
775  // Extract spectral model
776  m_spectral = xml_spectral(*spectral);
777 
778  // Optionally get temporal model
779  if (xml.elements("temporal")) {
780  const GXmlElement* temporal = xml.element("temporal", 0);
781  m_temporal = xml_temporal(*temporal);
782  }
783  else {
784  GModelTemporalConst constant;
785  m_temporal = constant.clone();
786  }
787 
788  // Read model attributes
789  read_attributes(xml);
790 
791  // Set parameter pointers
792  set_pointers();
793 
794  // Return
795  return;
796 }
797 
798 
799 /***********************************************************************//**
800  * @brief Write CTA effective area background model into XML element
801  *
802  * @param[in] xml XML element.
803  *
804  * Write CTA effective area background model information into an XML element.
805  * The XML element will have the following structure
806  *
807  * <source name="..." type="..." instrument="...">
808  * <spectrum type="...">
809  * ...
810  * </spectrum>
811  * </source>
812  *
813  * If the model contains a non-constant temporal model, the temporal
814  * component will also be written following the syntax
815  *
816  * <source name="..." type="..." instrument="...">
817  * <spectrum type="...">
818  * ...
819  * </spectrum>
820  * <temporal type="...">
821  * ...
822  * </temporal>
823  * </source>
824  *
825  * If no temporal component is found a constant model is assumed.
826  ***************************************************************************/
828 {
829  // Initialise pointer on source
830  GXmlElement* src = NULL;
831 
832  // Search corresponding source
833  int n = xml.elements("source");
834  for (int k = 0; k < n; ++k) {
835  GXmlElement* element = xml.element("source", k);
836  if (element->attribute("name") == name()) {
837  src = element;
838  break;
839  }
840  }
841 
842  // If we have a temporal model that is either not a constant, or a
843  // constant with a normalization value that differs from 1.0 then
844  // write the temporal component into the XML element. This logic
845  // assures compatibility with the Fermi/LAT format as this format
846  // does not handle temporal components.
847  bool write_temporal = ((m_temporal != NULL) &&
848  (m_temporal->type() != "Constant" ||
849  (*m_temporal)[0].value() != 1.0));
850 
851  // If no source with corresponding name was found then append one
852  if (src == NULL) {
853  src = xml.append("source");
854  if (spectral() != NULL) src->append(GXmlElement("spectrum"));
855  if (write_temporal) src->append(GXmlElement("temporal"));
856  }
857 
858  // Write spectral model
859  if (spectral() != NULL) {
860  GXmlElement* spec = src->element("spectrum", 0);
861  spectral()->write(*spec);
862  }
863 
864  // Optionally write temporal model
865  if (write_temporal) {
866  if (dynamic_cast<GModelTemporalConst*>(temporal()) == NULL) {
867  GXmlElement* temp = src->element("temporal", 0);
868  temporal()->write(*temp);
869  }
870  }
871 
872  // Write model attributes
873  write_attributes(*src);
874 
875  // Return
876  return;
877 }
878 
879 
880 /***********************************************************************//**
881  * @brief Print CTA effective area background model information
882  *
883  * @param[in] chatter Chattiness.
884  * @return String containing CTA effective area background model information.
885  ***************************************************************************/
886 std::string GCTAModelAeffBackground::print(const GChatter& chatter) const
887 {
888  // Initialise result string
889  std::string result;
890 
891  // Continue only if chatter is not silent
892  if (chatter != SILENT) {
893 
894  // Append header
895  result.append("=== GCTAModelAeffBackground ===");
896 
897  // Determine number of parameters per type
898  int n_spectral = (spectral() != NULL) ? spectral()->size() : 0;
899  int n_temporal = (temporal() != NULL) ? temporal()->size() : 0;
900 
901  // Append attributes
902  result.append("\n"+print_attributes());
903 
904  // Append model type
905  result.append("\n"+gammalib::parformat("Model type"));
906  if (n_spectral > 0) {
907  result.append("\""+spectral()->type()+"\"");
908  if (n_temporal > 0) {
909  result.append(" * ");
910  }
911  }
912  if (n_temporal > 0) {
913  result.append("\""+temporal()->type()+"\"");
914  }
915 
916  // Append parameters
917  result.append("\n"+gammalib::parformat("Number of parameters") +
918  gammalib::str(size()));
919  result.append("\n"+gammalib::parformat("Number of spectral par's") +
920  gammalib::str(n_spectral));
921  for (int i = 0; i < n_spectral; ++i) {
922  result.append("\n"+(*spectral())[i].print());
923  }
924  result.append("\n"+gammalib::parformat("Number of temporal par's") +
925  gammalib::str(n_temporal));
926  for (int i = 0; i < n_temporal; ++i) {
927  result.append("\n"+(*temporal())[i].print());
928  }
929 
930  } // endif: chatter was not silent
931 
932  // Return result
933  return result;
934 }
935 
936 
937 /*==========================================================================
938  = =
939  = Private methods =
940  = =
941  ==========================================================================*/
942 
943 /***********************************************************************//**
944  * @brief Initialise class members
945  ***************************************************************************/
947 {
948  // Initialise members
949  m_spectral = NULL;
950  m_temporal = NULL;
951  m_n_mc_energies = 100;
952 
953  // Initialise Npred cache
954  m_npred_names.clear();
955  m_npred_energies.clear();
956  m_npred_times.clear();
957  m_npred_values.clear();
958 
959  // Return
960  return;
961 }
962 
963 
964 /***********************************************************************//**
965  * @brief Copy class members
966  *
967  * @param[in] bgd CTA effective area background model.
968  ***************************************************************************/
970 {
971  // Copy cache
977 
978  // Clone spectral and temporal model components
979  m_spectral = (bgd.m_spectral != NULL) ? bgd.m_spectral->clone() : NULL;
980  m_temporal = (bgd.m_temporal != NULL) ? bgd.m_temporal->clone() : NULL;
981 
982  // Set parameter pointers
983  set_pointers();
984 
985  // Return
986  return;
987 }
988 
989 
990 /***********************************************************************//**
991  * @brief Delete class members
992  ***************************************************************************/
994 {
995  // Free memory
996  if (m_spectral != NULL) delete m_spectral;
997  if (m_temporal != NULL) delete m_temporal;
998 
999  // Signal free pointers
1000  m_spectral = NULL;
1001  m_temporal = NULL;
1002 
1003  // Return
1004  return;
1005 }
1006 
1007 
1008 /***********************************************************************//**
1009  * @brief Set pointers
1010  *
1011  * Set pointers to all model parameters. The pointers are stored in a vector
1012  * that is member of the GModelData base class.
1013  ***************************************************************************/
1015 {
1016  // Clear parameters
1017  m_pars.clear();
1018 
1019  // Determine the number of parameters
1020  int n_spectral = (spectral() != NULL) ? spectral()->size() : 0;
1021  int n_temporal = (temporal() != NULL) ? temporal()->size() : 0;
1022  int n_pars = n_spectral + n_temporal;
1023 
1024  // Continue only if there are parameters
1025  if (n_pars > 0) {
1026 
1027  // Gather spectral parameters
1028  for (int i = 0; i < n_spectral; ++i) {
1029  m_pars.push_back(&((*spectral())[i]));
1030  }
1031 
1032  // Gather temporal parameters
1033  for (int i = 0; i < n_temporal; ++i) {
1034  m_pars.push_back(&((*temporal())[i]));
1035  }
1036 
1037  }
1038 
1039  // Return
1040  return;
1041 }
1042 
1043 
1044 /***********************************************************************//**
1045  * @brief Verifies if model has all components
1046  *
1047  * Returns 'true' if models has a spectral and a temporal component.
1048  * Otherwise returns 'false'.
1049  ***************************************************************************/
1051 {
1052  // Set result
1053  bool result = ((spectral() != NULL) && (temporal() != NULL));
1054 
1055  // Return result
1056  return result;
1057 }
1058 
1059 
1060 /***********************************************************************//**
1061  * @brief Return pointer to spectral model from XML element
1062  *
1063  * @param[in] spectral XML element.
1064  * @return Pointer to spectral model.
1065  *
1066  * Returns pointer to spectral model that is defined in an XML element.
1067  ***************************************************************************/
1069 {
1070  // Get spectral model
1071  GModelSpectralRegistry registry;
1072  GModelSpectral* ptr = registry.alloc(spectral);
1073 
1074  // Return pointer
1075  return ptr;
1076 }
1077 
1078 
1079 /***********************************************************************//**
1080  * @brief Return pointer to temporal model from XML element
1081  *
1082  * @param[in] temporal XML element.
1083  * @return Pointer to temporal model.
1084  *
1085  * Returns pointer to temporal model that is defined in an XML element.
1086  ***************************************************************************/
1088 {
1089  // Get temporal model
1090  GModelTemporalRegistry registry;
1091  GModelTemporal* ptr = registry.alloc(temporal);
1092 
1093  // Return pointer
1094  return ptr;
1095 }
1096 
1097 
1098 /***********************************************************************//**
1099  * @brief Spatially integrate effective area for given energy
1100  *
1101  * @param[in] obs Observation.
1102  * @param[in] logE Log10 of reference energy in TeV.
1103  * @return Spatially integrated effective area for given energy.
1104  *
1105  * Spatially integrates the effective area for a given reference energy
1106  * over the region of interest.
1107  ***************************************************************************/
1109  const double& logE) const
1110 {
1111  // Set number of iterations for Romberg integration.
1112  static const int iter_theta = 6;
1113  static const int iter_phi = 6;
1114 
1115  // Get reference on CTA pointing, effective area response and event list
1116  // from observation
1117  const GCTAPointing& pnt = gammalib::cta_pnt(G_AEFF_INTEGRAL, obs);
1118  const GCTAAeff& aeff = gammalib::cta_rsp_aeff(G_AEFF_INTEGRAL, obs);
1120 
1121  // Get instrument direction of RoI centre
1122  GCTAInstDir roi_centre = pnt.instdir(events.roi().centre().dir());
1123 
1124  // Get ROI radius in radians
1125  double roi_radius = events.roi().radius() * gammalib::deg2rad;
1126 
1127  // Setup integration function
1129  logE,
1130  roi_centre,
1131  iter_phi);
1132 
1133  // Setup integration
1134  GIntegral integral(&integrand);
1135 
1136  // Set fixed number of iterations
1137  integral.fixed_iter(iter_theta);
1138 
1139  // Spatially integrate radial component
1140  double value = integral.romberg(0.0, roi_radius);
1141 
1142  // Debug: Check for NaN
1143  #if defined(G_NAN_CHECK)
1144  if (gammalib::is_notanumber(value) || gammalib::is_infinite(value)) {
1145  std::string origin = "GCTAModelAeffBackground::aeff_integral";
1146  std::string message = " NaN/Inf encountered (value=" +
1147  gammalib::str(value) + ", roi_radius=" +
1148  gammalib::str(roi_radius) + ")";
1149  gammalib::warning(origin, message);
1150  }
1151  #endif
1152 
1153  // Return
1154  return value;
1155 
1156 }
1157 
1158 
1159 /***********************************************************************//**
1160  * @brief Kernel for offset angle integration of effective area background model
1161  *
1162  * @param[in] theta Offset angle from ROI centre (radians).
1163  * @return Integration kernel value.
1164  *
1165  * Computes
1166  *
1167  * \f[
1168  * K(\rho | E, t) = \sin \theta \times
1169  * \int_{0}^{2\pi}
1170  * B(\theta,\phi | E, t) d\phi
1171  * \f]
1172  *
1173  * where \f$B(\theta,\phi | E, t)\f$ is the background model for a specific
1174  * observed energy \f$E\f$ and time \f$t\f$.
1175  ***************************************************************************/
1177 {
1178  // Initialise value
1179  double value = 0.0;
1180 
1181  // Continue only if offset angle is positive
1182  if (theta > 0.0) {
1183 
1184  // Setup phi integration kernel
1186  m_logE,
1187  m_roi_centre,
1188  theta);
1189 
1190  // Setup integration
1191  GIntegral integral(&integrand);
1192 
1193  // Set fixed number of iterations
1194  integral.fixed_iter(m_iter);
1195 
1196  // Integrate over phi
1197  value = integral.romberg(0.0, gammalib::twopi) * std::sin(theta);
1198 
1199  // Debug: Check for NaN
1200  #if defined(G_NAN_CHECK)
1201  if (gammalib::is_notanumber(value) || gammalib::is_infinite(value)) {
1202  std::string origin = "GCTAModelAeffBackground::npred_roi_kern_theta::eval"
1203  "(" + gammalib::str(theta) + ")";
1204  std::string message = " NaN/Inf encountered (value=" +
1205  gammalib::str(value) + ")";
1206  gammalib::warning(origin, message);
1207  }
1208  #endif
1209 
1210  } // endif: offset angle was positive
1211 
1212  // Return value
1213  return value;
1214 }
1215 
1216 
1217 /***********************************************************************//**
1218  * @brief Kernel for azimuth angle integration of effective area background
1219  * model
1220  *
1221  * @param[in] phi Azimuth angle around ROI centre (radians).
1222  * @return Integration kernel value.
1223  *
1224  * Computes
1225  *
1226  * \f[
1227  * B(\theta, \phi | E, t)
1228  * \f]
1229  ***************************************************************************/
1231 {
1232  // Compute detx and dety in radians
1233  double detx = m_roi_centre.detx();
1234  double dety = m_roi_centre.dety();
1235  if (m_theta > 0.0 ) {
1236  detx += m_theta * std::cos(phi);
1237  dety += m_theta * std::sin(phi);
1238  }
1239 
1240  // Convert into theta and phi
1241  GCTAInstDir dir(detx, dety);
1242  double theta_prime = dir.theta();
1243  double phi_prime = dir.phi();
1244 
1245  // Get background value
1246  double value = (*m_aeff)(m_logE, theta_prime, phi_prime);
1247 
1248  // Debug: Check for NaN
1249  #if defined(G_NAN_CHECK)
1250  if (gammalib::is_notanumber(value) || gammalib::is_infinite(value)) {
1251  std::string origin = "GCTAModelAeffBackground::npred_roi_kern_phi::eval"
1252  "(" + gammalib::str(phi) + ")";
1253  std::string message = " NaN/Inf encountered (value=" +
1254  gammalib::str(value) + ", theta=" +
1255  gammalib::str(m_theta) + ")";
1256  gammalib::warning(origin, message);
1257  }
1258  #endif
1259 
1260  // Return Npred
1261  return value;
1262 }
virtual void write(GXmlElement &xml) const =0
CTA background model base class definition.
Constant temporal model class interface definition.
void init_members(void)
Initialise class members.
#define G_MC
CTA Aeff background model class definition.
double eval(const double &phi)
Kernel for azimuth angle integration of effective area background model.
virtual GTimes mc(const double &rate, const GTime &tmin, const GTime &tmax, GRan &ran) const =0
double norm(const GVector &vector)
Computes vector norm.
Definition: GVector.cpp:864
double romberg(std::vector< double > bounds, const int &order=5)
Perform Romberg integration.
Definition: GIntegral.cpp:382
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:770
void warning(const std::string &origin, const std::string &message)
Emits warning.
Definition: GTools.cpp:1386
Abstract spectral model base class.
Interface definition for the model registry class.
virtual std::string instrument(void) const =0
Abstract temporal model base class.
virtual void roi(const GRoi &roi)
Set ROI.
const GCTAModelAeffBackground g_cta_aeff_background_seed
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:163
#define G_AEFF_INTEGRAL
GVector cos(const GVector &vector)
Computes cosine of vector elements.
Definition: GVector.cpp:1190
virtual double eval(const GEvent &event, const GObservation &obs, const bool &gradients=false) const
Return background rate in units of events MeV s sr .
virtual void read(const GXmlElement &xml)
Read CTA effective area background model from XML element.
virtual double eval(const GEnergy &srcEng, const GTime &srcTime=GTime(), const bool &gradients=false) const =0
GModelTemporal * xml_temporal(const GXmlElement &temporal) const
Return pointer to temporal model from XML element.
void append(const GCTAEventAtom &event)
Append event to event list.
GModelTemporal * temporal(void) const
Return temporal model component.
CTA instrument response function class definition.
const GCTAPointing & cta_pnt(const std::string &origin, const GObservation &obs)
Retrieve CTA pointing from generic observation.
Abstract interface for the event classes.
Definition: GEvent.hpp:71
CTA event list class.
const double & radius(void) const
Returns radius of region of interest in degrees.
Definition: GCTARoi.hpp:125
XML element node class.
Definition: GXmlElement.hpp:48
void gti(const GGti &gti)
Set Good Time Intervals.
Definition: GEvents.cpp:154
const GCTAAeff * m_aeff
Pointer to effectve area.
void free_members(void)
Delete class members.
Definition: GModelData.cpp:300
Spectral model registry class definition.
const double & zenith(void) const
Return pointing zenith angle.
virtual void write(GXmlElement &xml) const =0
double acos(const double &arg)
Computes acos by avoiding NaN due to rounding errors.
Definition: GMath.cpp:69
double log10TeV(void) const
Return log10 of energy in TeV.
Definition: GEnergy.cpp:457
GModelSpectral * m_spectral
Spectral model.
Random number generator class.
Definition: GRan.hpp:44
virtual std::string type(void) const
Return data model type.
Spectral nodes model class definition.
Interface for the CTA region of interest class.
Definition: GCTARoi.hpp:49
virtual int elements(void) const
Return number of GXMLElement children of node.
Definition: GXmlNode.cpp:586
Time class.
Definition: GTime.hpp:55
Gammalib tools definition.
GIntegral class interface definition.
Definition: GIntegral.hpp:46
virtual GModelTemporalConst * clone(void) const
Clone constant temporal model.
virtual const GTime & time(void) const =0
void dir(const GSkyDir &dir)
Set sky direction.
double eval(const double &theta)
Kernel for offset angle integration of effective area background model.
int size(void) const
Return number of Good Time Intervals.
Definition: GGti.hpp:154
Spectral nodes model class.
GCTAModelAeffBackground(void)
Void constructor.
void id(const std::string &id)
Set observation identifier.
Time container class.
Definition: GTimes.hpp:45
bool is_notanumber(const double &x)
Signal if argument is not a number.
Definition: GTools.hpp:201
std::vector< GModelPar * > m_pars
Pointers to all model parameters.
Definition: GModel.hpp:178
int size(void) const
Return number of parameters in model.
Definition: GModel.hpp:233
bool is_infinite(const double &x)
Signal if argument is infinite.
Definition: GTools.hpp:184
Definition of support function used by CTA classes.
virtual GCTAModelAeffBackground * clone(void) const
Clone CTA effective area background model.
const double deg2rad
Definition: GMath.hpp:43
virtual GCTAEventList * mc(const GObservation &obs, GRan &ran) const
Return simulated list of events.
std::vector< double > m_npred_values
Model values.
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:261
void reserve(const int &number)
Reserves space for events.
virtual void clear(void)
Clear CTA effective area background model.
double theta(void) const
Return offset angle (in radians)
virtual double npred(const GEnergy &obsEng, const GTime &obsTime, const GObservation &obs) const
Return spatially integrated background rate in units of events MeV s .
void rotate_deg(const double &phi, const double &theta)
Rotate sky direction by zenith and azimuth angle.
Definition: GSkyDir.cpp:424
const GEnergy & emin(void) const
Return minimum energy of all intervals.
Definition: GEbounds.hpp:187
std::vector< GEnergy > m_npred_energies
Model energy.
const GCTAInstDir & centre(void) const
Returns region of interest centre.
Definition: GCTARoi.hpp:111
int size(void) const
Return number of parameters.
GCTAInstDir instdir(const GSkyDir &skydir) const
Get instrument direction from sky direction.
GModelTemporal * alloc(const GXmlElement &xml) const
Allocate temporal model that is found in XML element.
virtual const GEnergy & energy(void) const =0
double uniform(void)
Returns random double precision floating value in range 0 to 1.
Definition: GRan.cpp:242
CTA pointing class.
Abstract data model class.
Definition: GModelData.hpp:55
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 std::string print(const GChatter &chatter=NORMAL) const
Print CTA effective area background model information.
virtual GModelData & operator=(const GModelData &model)
Assignment operator.
Definition: GModelData.cpp:125
GChatter
Definition: GTypemaps.hpp:33
std::vector< GTime > m_npred_times
Model time.
std::vector< std::string > m_npred_names
Model names.
Good Time Interval class.
Definition: GGti.hpp:62
void copy_members(const GCTAModelAeffBackground &bgd)
Copy class members.
virtual double eval(const GTime &srcTime, const bool &gradients=false) const =0
Abstract observation base class.
virtual double max(const double &logE, const double &zenith, const double &azimuth, const bool &etrue=true) const =0
const GTime & tstop(void) const
Returns latest stop time in Good Time Intervals.
Definition: GGti.hpp:207
const GTime & tstart(void) const
Returns earliest start time in Good Time Intervals.
Definition: GGti.hpp:194
CTA observation class 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:723
virtual GModelSpectral * clone(void) const =0
Clones object.
Abstract base class for the CTA effective area.
Definition: GCTAAeff.hpp:54
const double & azimuth(void) const
Return pointing azimuth angle.
virtual ~GCTAModelAeffBackground(void)
Destructor.
virtual std::string type(void) const =0
int iter_phi(const double &rho, const double &resolution, const int &iter_min, const int &iter_max)
Determine number of azimuthal Romberg iterations.
virtual GXmlElement * element(const int &index)
Return pointer to GXMLElement child.
Definition: GXmlNode.cpp:640
Temporal model registry class definition.
GModelTemporal * m_temporal
Temporal model.
GModelSpectral * spectral(void) const
Return spectral model component.
void free_members(void)
Delete class members.
bool valid_model(void) const
Verifies if model has all components.
CTA event atom class.
GVector sin(const GVector &vector)
Computes sine of vector elements.
Definition: GVector.cpp:1316
double aeff_integral(const GObservation &obs, const double &logE) const
Spatially integrate effective area for given energy.
GCTAModelAeffBackground & operator=(const GCTAModelAeffBackground &bgd)
Assignment operator.
double dist(const GSkyDir &dir) const
Compute angular distance between sky directions in radians.
Definition: GSkyDir.hpp:271
void fixed_iter(const int &iter)
Set fixed number of iterations.
Definition: GIntegral.hpp:189
virtual void write(GXmlElement &xml) const
Write CTA effective area background model into XML element.
virtual GModelTemporal * clone(void) const =0
Clones object.
void read_attributes(const GXmlElement &xml)
Read model attributes.
Definition: GModel.cpp:684
GModelSpectral * xml_spectral(const GXmlElement &spectral) const
Return pointer to spectral model from XML element.
CTA instrument direction class.
Definition: GCTAInstDir.hpp:63
const GSkyDir & dir(void) const
Return pointing sky direction.
const double twopi
Definition: GMath.hpp:36
const GEnergy & emax(void) const
Return maximum energy of all intervals.
Definition: GEbounds.hpp:199
Model registry class definition.
const double rad2deg
Definition: GMath.hpp:44
virtual GXmlNode * append(const GXmlNode &node)
Append XML child node.
Definition: GXmlNode.cpp:287
std::string parformat(const std::string &s, const int &indent=0)
Convert string in parameter format.
Definition: GTools.cpp:1143
const GCTAAeff & cta_rsp_aeff(const std::string &origin, const GObservation &obs)
Retrieve CTA effective area response from generic observation.
const GTime & tstart(void) const
Return start time.
Definition: GEvents.hpp:146
double phi(void) const
Return azimuth angle (in radians)
Integration class interface definition.
int m_n_mc_energies
Energy sampling for MC spectrum.
Sky direction class.
Definition: GSkyDir.hpp:62
#define G_EVAL
void init_members(void)
Initialise class members.
Definition: GModelData.cpp:278
Constant temporal model class.
const GCTAInstDir & dir(void) const
Return instrument direction.
Mathematical function definitions.
Class that handles energies in a unit independent way.
Definition: GEnergy.hpp:48
std::string str(const unsigned short int &value)
Convert unsigned short integer value into string.
Definition: GTools.cpp:489
void set_pointers(void)
Set pointers.