GammaLib  1.7.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 =
567  GEbounds(m_n_mc_energies, ebounds.emin(), ebounds.emax(), true);
569  for (int i = 0; i < spectral_ebounds.size(); ++i) {
570  GEnergy energy = spectral_ebounds.elogmean(i);
571  double intensity = aeff_integral(obs, energy.log10TeV());
572  double norm = m_spectral->eval(energy, events.tstart());
573  double arg = norm * intensity;
574  if (arg > 0.0) {
575  spectral.append(energy, arg);
576  }
577  }
578 
579  // Loop over all energy bins
580  for (int ieng = 0; ieng < ebounds.size(); ++ieng) {
581 
582  // Compute the background rate in model within the energy
583  // boundaries from spectral component (units: cts/s).
584  // Note that the time here is ontime. Deadtime correction will
585  // be done later.
586  double rate = spectral.flux(ebounds.emin(ieng), ebounds.emax(ieng));
587 
588  // Debug option: dump rate
589  #if defined(G_DUMP_MC)
590  std::cout << "GCTAModelAeffBackground::mc(\"" << name() << "\": ";
591  std::cout << "rate=" << rate << " cts/s)" << std::endl;
592  #endif
593 
594  // If the rate is not positive then skip this energy bins
595  if (rate <= 0.0) {
596  continue;
597  }
598 
599  // Loop over all good time intervals
600  for (int itime = 0; itime < gti.size(); ++itime) {
601 
602  // Get Monte Carlo event arrival times from temporal model
603  GTimes times = m_temporal->mc(rate,
604  gti.tstart(itime),
605  gti.tstop(itime),
606  ran);
607 
608  // Get number of events
609  int n_events = times.size();
610 
611  // Reserve space for events
612  if (n_events > 0) {
613  list->reserve(n_events);
614  }
615 
616  // Debug option: provide number of times and initialize
617  // statisics
618  #if defined(G_DUMP_MC)
619  std::cout << " Interval " << itime;
620  std::cout << " times=" << n_events << std::endl;
621  int n_killed_by_roi = 0;
622  #endif
623 
624  // Loop over events
625  for (int i = 0; i < n_events; ++i) {
626 
627  // Get Monte Carlo event energy from spectral model
628  GEnergy energy = spectral.mc(ebounds.emin(ieng),
629  ebounds.emax(ieng),
630  times[i],
631  ran);
632 
633  // Get maximum effective area for rejection method
634  double max_aeff = aeff.max(energy.log10TeV(), pnt.zenith(),
635  pnt.azimuth(), false);
636 
637  // Skip event if the maximum effective area is not positive
638  if (max_aeff <= 0.0) {
639  continue;
640  }
641 
642  // Initialise randomised coordinates
643  double offset = 0.0;
644  double phi = 0.0;
645 
646  // Initialise acceptance fraction and counter of zeros for
647  // rejection method
648  double acceptance_fraction = 0.0;
649  int zeros = 0;
650 
651  // Start rejection method loop
652  do {
653 
654  // Throw random offset and azimuth angle in
655  // considered range
656  offset = std::acos(1.0 - ran.uniform() *
657  (1.0 - cos_max_theta));
658  phi = ran.uniform() * gammalib::twopi;
659 
660  // Compute function value at this offset angle
661  double value = aeff(energy.log10TeV(), offset, phi,
662  pnt.zenith(), pnt.azimuth(),
663  false);
664 
665  // If the value is not positive then increment the
666  // zeros counter and fall through. The counter assures
667  // that this loop does not lock up.
668  if (value <= 0.0) {
669  zeros++;
670  continue;
671  }
672 
673  // Value is non-zero so reset the zeros counter
674  zeros = 0;
675 
676  // Compute acceptance fraction
677  acceptance_fraction = value / max_aeff;
678 
679  } while ((ran.uniform() > acceptance_fraction) &&
680  (zeros < 1000));
681 
682  // If the zeros counter is non-zero then the loop was
683  // exited due to exhaustion and the event is skipped
684  if (zeros > 0) {
685  continue;
686  }
687 
688  // Rotate pointing direction by offset and azimuth angle
689  GSkyDir skydir = pnt.dir();
690  skydir.rotate_deg(phi * gammalib::rad2deg,
691  offset * gammalib::rad2deg);
692 
693  // Convert rotated pointing direction in instrument system
694  GCTAInstDir mc_dir = pnt.instdir(skydir);
695 
696  // Allocate event
697  GCTAEventAtom event;
698 
699  // Set event attributes
700  event.dir(mc_dir);
701  event.energy(energy);
702  event.time(times[i]);
703 
704  // Append event to list if it falls in ROI
705  if (events.roi().contains(event)) {
706  list->append(event);
707  }
708  #if defined(G_DUMP_MC)
709  else {
710  n_killed_by_roi++;
711  }
712  #endif
713 
714  } // endfor: looped over all events
715 
716  // Debug option: provide statisics
717  #if defined(G_DUMP_MC)
718  std::cout << " Killed by ROI=";
719  std::cout << n_killed_by_roi << std::endl;
720  #endif
721 
722  } // endfor: looped over all GTIs
723 
724  } // endfor: looped over all energy boundaries
725 
726  } // endif: model was valid
727 
728  // Return
729  return list;
730 }
731 
732 
733 /***********************************************************************//**
734  * @brief Read CTA effective area background model from XML element
735  *
736  * @param[in] xml XML element.
737  *
738  * Set up CTA effective area background model from the information provided by
739  * an XML element. The XML element is expected to have the following
740  * structure
741  *
742  * <source name="..." type="..." instrument="...">
743  * <spectrum type="...">
744  * ...
745  * </spectrum>
746  * </source>
747  *
748  * Optionally, a temporal model may be provided using the following
749  * syntax
750  *
751  * <source name="..." type="..." instrument="...">
752  * <spectrum type="...">
753  * ...
754  * </spectrum>
755  * <temporalModel type="...">
756  * ...
757  * </temporalModel>
758  * </source>
759  *
760  * If no temporal component is found a constant model is assumed.
761  ***************************************************************************/
763 {
764  // Clear model
765  clear();
766 
767  // Initialise XML elements
768  const GXmlElement* spectral = NULL;
769  const GXmlElement* temporal = NULL;
770 
771  // Get pointer on spectrum
772  spectral = xml.element("spectrum", 0);
773 
774  // Extract spectral model
775  m_spectral = xml_spectral(*spectral);
776 
777  // Optionally get temporal model
778  if (xml.elements("temporalModel")) {
779  temporal = xml.element("temporalModel", 0);
780  m_temporal = xml_temporal(*temporal);
781  }
782  else {
783  GModelTemporalConst constant;
784  m_temporal = constant.clone();
785  }
786 
787  // Read model attributes
788  read_attributes(xml);
789 
790  // Set parameter pointers
791  set_pointers();
792 
793  // Return
794  return;
795 }
796 
797 
798 /***********************************************************************//**
799  * @brief Write CTA effective area background model into XML element
800  *
801  * @param[in] xml XML element.
802  *
803  * Write CTA effective area background model information into an XML element.
804  * The XML element will have the following structure
805  *
806  * <source name="..." type="..." instrument="...">
807  * <spectrum type="...">
808  * ...
809  * </spectrum>
810  * </source>
811  *
812  * If the model contains a non-constant temporal model, the temporal
813  * component will also be written following the syntax
814  *
815  * <source name="..." type="..." instrument="...">
816  * <spectrum type="...">
817  * ...
818  * </spectrum>
819  * <temporalModel type="...">
820  * ...
821  * </temporalModel>
822  * </source>
823  *
824  * If no temporal component is found a constant model is assumed.
825  ***************************************************************************/
827 {
828  // Initialise pointer on source
829  GXmlElement* src = NULL;
830 
831  // Search corresponding source
832  int n = xml.elements("source");
833  for (int k = 0; k < n; ++k) {
834  GXmlElement* element = xml.element("source", k);
835  if (element->attribute("name") == name()) {
836  src = element;
837  break;
838  }
839  }
840 
841  // If we have a temporal model that is either not a constant, or a
842  // constant with a normalization value that differs from 1.0 then
843  // write the temporal component into the XML element. This logic
844  // assures compatibility with the Fermi/LAT format as this format
845  // does not handle temporal components.
846  bool write_temporal = ((m_temporal != NULL) &&
847  (m_temporal->type() != "Constant" ||
848  (*m_temporal)[0].value() != 1.0));
849 
850  // If no source with corresponding name was found then append one
851  if (src == NULL) {
852  src = xml.append("source");
853  if (spectral() != NULL) src->append(GXmlElement("spectrum"));
854  if (write_temporal) src->append(GXmlElement("temporalModel"));
855  }
856 
857  // Write spectral model
858  if (spectral() != NULL) {
859  GXmlElement* spec = src->element("spectrum", 0);
860  spectral()->write(*spec);
861  }
862 
863  // Optionally write temporal model
864  if (write_temporal) {
865  if (dynamic_cast<GModelTemporalConst*>(temporal()) == NULL) {
866  GXmlElement* temp = src->element("temporalModel", 0);
867  temporal()->write(*temp);
868  }
869  }
870 
871  // Write model attributes
872  write_attributes(*src);
873 
874  // Return
875  return;
876 }
877 
878 
879 /***********************************************************************//**
880  * @brief Print CTA effective area background model information
881  *
882  * @param[in] chatter Chattiness.
883  * @return String containing CTA effective area background model information.
884  ***************************************************************************/
885 std::string GCTAModelAeffBackground::print(const GChatter& chatter) const
886 {
887  // Initialise result string
888  std::string result;
889 
890  // Continue only if chatter is not silent
891  if (chatter != SILENT) {
892 
893  // Append header
894  result.append("=== GCTAModelAeffBackground ===");
895 
896  // Determine number of parameters per type
897  int n_spectral = (spectral() != NULL) ? spectral()->size() : 0;
898  int n_temporal = (temporal() != NULL) ? temporal()->size() : 0;
899 
900  // Append attributes
901  result.append("\n"+print_attributes());
902 
903  // Append model type
904  result.append("\n"+gammalib::parformat("Model type"));
905  if (n_spectral > 0) {
906  result.append("\""+spectral()->type()+"\"");
907  if (n_temporal > 0) {
908  result.append(" * ");
909  }
910  }
911  if (n_temporal > 0) {
912  result.append("\""+temporal()->type()+"\"");
913  }
914 
915  // Append parameters
916  result.append("\n"+gammalib::parformat("Number of parameters") +
917  gammalib::str(size()));
918  result.append("\n"+gammalib::parformat("Number of spectral par's") +
919  gammalib::str(n_spectral));
920  for (int i = 0; i < n_spectral; ++i) {
921  result.append("\n"+(*spectral())[i].print());
922  }
923  result.append("\n"+gammalib::parformat("Number of temporal par's") +
924  gammalib::str(n_temporal));
925  for (int i = 0; i < n_temporal; ++i) {
926  result.append("\n"+(*temporal())[i].print());
927  }
928 
929  } // endif: chatter was not silent
930 
931  // Return result
932  return result;
933 }
934 
935 
936 /*==========================================================================
937  = =
938  = Private methods =
939  = =
940  ==========================================================================*/
941 
942 /***********************************************************************//**
943  * @brief Initialise class members
944  ***************************************************************************/
946 {
947  // Initialise members
948  m_spectral = NULL;
949  m_temporal = NULL;
950  m_n_mc_energies = 100;
951 
952  // Initialise Npred cache
953  m_npred_names.clear();
954  m_npred_energies.clear();
955  m_npred_times.clear();
956  m_npred_values.clear();
957 
958  // Return
959  return;
960 }
961 
962 
963 /***********************************************************************//**
964  * @brief Copy class members
965  *
966  * @param[in] bgd CTA effective area background model.
967  ***************************************************************************/
969 {
970  // Copy cache
976 
977  // Clone spectral and temporal model components
978  m_spectral = (bgd.m_spectral != NULL) ? bgd.m_spectral->clone() : NULL;
979  m_temporal = (bgd.m_temporal != NULL) ? bgd.m_temporal->clone() : NULL;
980 
981  // Set parameter pointers
982  set_pointers();
983 
984  // Return
985  return;
986 }
987 
988 
989 /***********************************************************************//**
990  * @brief Delete class members
991  ***************************************************************************/
993 {
994  // Free memory
995  if (m_spectral != NULL) delete m_spectral;
996  if (m_temporal != NULL) delete m_temporal;
997 
998  // Signal free pointers
999  m_spectral = NULL;
1000  m_temporal = NULL;
1001 
1002  // Return
1003  return;
1004 }
1005 
1006 
1007 /***********************************************************************//**
1008  * @brief Set pointers
1009  *
1010  * Set pointers to all model parameters. The pointers are stored in a vector
1011  * that is member of the GModelData base class.
1012  ***************************************************************************/
1014 {
1015  // Clear parameters
1016  m_pars.clear();
1017 
1018  // Determine the number of parameters
1019  int n_spectral = (spectral() != NULL) ? spectral()->size() : 0;
1020  int n_temporal = (temporal() != NULL) ? temporal()->size() : 0;
1021  int n_pars = n_spectral + n_temporal;
1022 
1023  // Continue only if there are parameters
1024  if (n_pars > 0) {
1025 
1026  // Gather spectral parameters
1027  for (int i = 0; i < n_spectral; ++i) {
1028  m_pars.push_back(&((*spectral())[i]));
1029  }
1030 
1031  // Gather temporal parameters
1032  for (int i = 0; i < n_temporal; ++i) {
1033  m_pars.push_back(&((*temporal())[i]));
1034  }
1035 
1036  }
1037 
1038  // Return
1039  return;
1040 }
1041 
1042 
1043 /***********************************************************************//**
1044  * @brief Verifies if model has all components
1045  *
1046  * Returns 'true' if models has a spectral and a temporal component.
1047  * Otherwise returns 'false'.
1048  ***************************************************************************/
1050 {
1051  // Set result
1052  bool result = ((spectral() != NULL) && (temporal() != NULL));
1053 
1054  // Return result
1055  return result;
1056 }
1057 
1058 
1059 /***********************************************************************//**
1060  * @brief Return pointer to spectral model from XML element
1061  *
1062  * @param[in] spectral XML element.
1063  * @return Pointer to spectral model.
1064  *
1065  * Returns pointer to spectral model that is defined in an XML element.
1066  ***************************************************************************/
1068 {
1069  // Get spectral model
1070  GModelSpectralRegistry registry;
1071  GModelSpectral* ptr = registry.alloc(spectral);
1072 
1073  // Return pointer
1074  return ptr;
1075 }
1076 
1077 
1078 /***********************************************************************//**
1079  * @brief Return pointer to temporal model from XML element
1080  *
1081  * @param[in] temporal XML element.
1082  * @return Pointer to temporal model.
1083  *
1084  * Returns pointer to temporal model that is defined in an XML element.
1085  ***************************************************************************/
1087 {
1088  // Get temporal model
1089  GModelTemporalRegistry registry;
1090  GModelTemporal* ptr = registry.alloc(temporal);
1091 
1092  // Return pointer
1093  return ptr;
1094 }
1095 
1096 
1097 /***********************************************************************//**
1098  * @brief Spatially integrate effective area for given energy
1099  *
1100  * @param[in] obs Observation.
1101  * @param[in] logE Log10 of reference energy in TeV.
1102  * @return Spatially integrated effective area for given energy.
1103  *
1104  * Spatially integrates the effective area for a given reference energy
1105  * over the region of interest.
1106  ***************************************************************************/
1108  const double& logE) const
1109 {
1110  // Initialise result
1111  double value = 0.0;
1112 
1113  // Set number of iterations for Romberg integration.
1114  static const int iter_theta = 6;
1115  static const int iter_phi = 6;
1116 
1117  // Get reference on CTA pointing, effective area response and event list
1118  // from observation
1119  const GCTAPointing& pnt = gammalib::cta_pnt(G_AEFF_INTEGRAL, obs);
1120  const GCTAAeff& aeff = gammalib::cta_rsp_aeff(G_AEFF_INTEGRAL, obs);
1122 
1123  // Get instrument direction of RoI centre
1124  GCTAInstDir roi_centre = pnt.instdir(events.roi().centre().dir());
1125 
1126  // Get ROI radius in radians
1127  double roi_radius = events.roi().radius() * gammalib::deg2rad;
1128 
1129  // Setup integration function
1131  logE,
1132  roi_centre,
1133  iter_phi);
1134 
1135  // Setup integration
1136  GIntegral integral(&integrand);
1137 
1138  // Set fixed number of iterations
1139  integral.fixed_iter(iter_theta);
1140 
1141  // Spatially integrate radial component
1142  value = integral.romberg(0.0, roi_radius);
1143 
1144  // Debug: Check for NaN
1145  #if defined(G_NAN_CHECK)
1146  if (gammalib::is_notanumber(value) || gammalib::is_infinite(value)) {
1147  std::string origin = "GCTAModelAeffBackground::aeff_integral";
1148  std::string message = " NaN/Inf encountered (value=" +
1149  gammalib::str(value) + ", roi_radius=" +
1150  gammalib::str(roi_radius) + ")";
1151  gammalib::warning(origin, message);
1152  }
1153  #endif
1154 
1155  // Return
1156  return value;
1157 
1158 }
1159 
1160 
1161 /***********************************************************************//**
1162  * @brief Kernel for offset angle integration of effective area background model
1163  *
1164  * @param[in] theta Offset angle from ROI centre (radians).
1165  * @return Integration kernel value.
1166  *
1167  * Computes
1168  *
1169  * \f[
1170  * K(\rho | E, t) = \sin \theta \times
1171  * \int_{0}^{2\pi}
1172  * B(\theta,\phi | E, t) d\phi
1173  * \f]
1174  *
1175  * where \f$B(\theta,\phi | E, t)\f$ is the background model for a specific
1176  * observed energy \f$E\f$ and time \f$t\f$.
1177  ***************************************************************************/
1179 {
1180  // Initialise value
1181  double value = 0.0;
1182 
1183  // Continue only if offset angle is positive
1184  if (theta > 0.0) {
1185 
1186  // Setup phi integration kernel
1188  m_logE,
1189  m_roi_centre,
1190  theta);
1191 
1192  // Setup integration
1193  GIntegral integral(&integrand);
1194 
1195  // Set fixed number of iterations
1196  integral.fixed_iter(m_iter);
1197 
1198  // Integrate over phi
1199  value = integral.romberg(0.0, gammalib::twopi) * std::sin(theta);
1200 
1201  // Debug: Check for NaN
1202  #if defined(G_NAN_CHECK)
1203  if (gammalib::is_notanumber(value) || gammalib::is_infinite(value)) {
1204  std::string origin = "GCTAModelAeffBackground::npred_roi_kern_theta::eval"
1205  "(" + gammalib::str(theta) + ")";
1206  std::string message = " NaN/Inf encountered (value=" +
1207  gammalib::str(value) + ")";
1208  gammalib::warning(origin, message);
1209  }
1210  #endif
1211 
1212  } // endif: offset angle was positive
1213 
1214  // Return value
1215  return value;
1216 }
1217 
1218 
1219 /***********************************************************************//**
1220  * @brief Kernel for azimuth angle integration of effective area background
1221  * model
1222  *
1223  * @param[in] phi Azimuth angle around ROI centre (radians).
1224  * @return Integration kernel value.
1225  *
1226  * Computes
1227  *
1228  * \f[
1229  * B(\theta, \phi | E, t)
1230  * \f]
1231  ***************************************************************************/
1233 {
1234  // Compute detx and dety in radians
1235  double detx = m_roi_centre.detx();
1236  double dety = m_roi_centre.dety();
1237  if (m_theta > 0.0 ) {
1238  detx += m_theta * std::cos(phi);
1239  dety += m_theta * std::sin(phi);
1240  }
1241 
1242  // Convert into theta and phi
1243  GCTAInstDir dir(detx, dety);
1244  double theta_prime = dir.theta();
1245  double phi_prime = dir.phi();
1246 
1247  // Get background value
1248  double value = (*m_aeff)(m_logE, theta_prime, phi_prime);
1249 
1250  // Debug: Check for NaN
1251  #if defined(G_NAN_CHECK)
1252  if (gammalib::is_notanumber(value) || gammalib::is_infinite(value)) {
1253  std::string origin = "GCTAModelAeffBackground::npred_roi_kern_phi::eval"
1254  "(" + gammalib::str(phi) + ")";
1255  std::string message = " NaN/Inf encountered (value=" +
1256  gammalib::str(value) + ", theta=" +
1257  gammalib::str(m_theta) + ")";
1258  gammalib::warning(origin, message);
1259  }
1260  #endif
1261 
1262  // Return Npred
1263  return value;
1264 }
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:821
double romberg(std::vector< double > bounds, const int &order=5)
Perform Romberg integration.
Definition: GIntegral.cpp:380
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:748
void warning(const std::string &origin, const std::string &message)
Emits warning.
Definition: GTools.cpp:1265
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:157
#define G_AEFF_INTEGRAL
GVector cos(const GVector &vector)
Computes cosine of vector elements.
Definition: GVector.cpp:1100
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:124
XML element node class.
Definition: GXmlElement.hpp:47
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:181
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: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.
GIntegral class interface definition.
Definition: GIntegral.hpp:45
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:153
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:184
std::vector< GModelPar * > m_pars
Pointers to all model parameters.
Definition: GModel.hpp:169
int size(void) const
Return number of parameters in model.
Definition: GModel.hpp:217
bool is_infinite(const double &x)
Signal if argument is infinite.
Definition: GTools.hpp:167
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:245
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:321
const GEnergy & emin(void) const
Return minimum energy of all intervals.
Definition: GEbounds.hpp:181
std::vector< GEnergy > m_npred_energies
Model energy.
const GCTAInstDir & centre(void) const
Returns region of interest centre.
Definition: GCTARoi.hpp:110
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:53
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:120
GChatter
Definition: GTypemaps.hpp:33
std::vector< GTime > m_npred_times
Model time.
std::vector< std::string > m_npred_names
Model names.
GEnergy elogmean(const int &index) const
Returns logarithmic mean energy for a given energy interval.
Definition: GEbounds.cpp:1140
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:206
const GTime & tstart(void) const
Returns earliest start time in Good Time Intervals.
Definition: GGti.hpp:193
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:704
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:634
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:1226
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:265
void fixed_iter(const int &iter)
Set fixed number of iterations.
Definition: GIntegral.hpp:181
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:668
GModelSpectral * xml_spectral(const GXmlElement &spectral) const
Return pointer to spectral model from XML element.
CTA instrument direction class.
Definition: GCTAInstDir.hpp:59
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:193
Model registry class definition.
const double rad2deg
Definition: GMath.hpp:44
virtual GXmlNode * append(const GXmlNode &node)
Append XML child node.
Definition: GXmlNode.cpp:284
std::string parformat(const std::string &s, const int &indent=0)
Convert string in parameter format.
Definition: GTools.cpp:1022
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:159
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:413
void set_pointers(void)
Set pointers.