GammaLib  1.7.0.dev
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
GCTAModelIrfBackground.cpp
Go to the documentation of this file.
1 /***************************************************************************
2  * GCTAModelIrfBackground.cpp - CTA IRF background model class *
3  * ----------------------------------------------------------------------- *
4  * copyright (C) 2014-2020 by Juergen Knoedlseder *
5  * ----------------------------------------------------------------------- *
6  * *
7  * This program is free software: you can redistribute it and/or modify *
8  * it under the terms of the GNU General Public License as published by *
9  * the Free Software Foundation, either version 3 of the License, or *
10  * (at your option) any later version. *
11  * *
12  * This program is distributed in the hope that it will be useful, *
13  * but WITHOUT ANY WARRANTY; without even the implied warranty of *
14  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
15  * GNU General Public License for more details. *
16  * *
17  * You should have received a copy of the GNU General Public License *
18  * along with this program. If not, see <http://www.gnu.org/licenses/>. *
19  * *
20  ***************************************************************************/
21 /**
22  * @file GCTAModelIrfBackground.cpp
23  * @brief CTA IRF background model class implementation
24  * @author Juergen Knoedlseder
25  */
26 
27 /* __ Includes ___________________________________________________________ */
28 #ifdef HAVE_CONFIG_H
29 #include <config.h>
30 #endif
31 #include "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_inst_background_registry(&g_cta_inst_background_seed);
48 
49 /* __ Method name definitions ____________________________________________ */
50 #define G_EVAL "GCTAModelIrfBackground::eval(GEvent&, GObservation&, bool&)"
51 #define G_NPRED "GCTAModelIrfBackground::npred(GEnergy&, GTime&,"\
52  " GObservation&)"
53 #define G_MC "GCTAModelIrfBackground::mc(GObservation&, GRan&)"
54 #define G_XML_SPECTRAL "GCTAModelIrfBackground::xml_spectral(GXmlElement&)"
55 #define G_XML_TEMPORAL "GCTAModelIrfBackground::xml_temporal(GXmlElement&)"
56 
57 /* __ Macros _____________________________________________________________ */
58 
59 /* __ Coding definitions _________________________________________________ */
60 #define G_USE_NPRED_CACHE
61 
62 /* __ Debug definitions __________________________________________________ */
63 //#define G_DUMP_MC
64 //#define G_DEBUG_NPRED
65 
66 /* __ Constants __________________________________________________________ */
67 
68 
69 /*==========================================================================
70  = =
71  = Constructors/destructors =
72  = =
73  ==========================================================================*/
74 
75 /***********************************************************************//**
76  * @brief Void constructor
77  ***************************************************************************/
79 {
80  // Initialise class members
81  init_members();
82 
83  // Return
84  return;
85 }
86 
87 
88 /***********************************************************************//**
89  * @brief Constructor
90  *
91  * @param[in] xml XML element.
92  *
93  * Constructs a CTA instrumental background model from the information
94  * provided by an XML element. The XML element is expected to have the
95  * following structure
96  *
97  * <source name="..." type="..." instrument="...">
98  * <spectrum type="...">
99  * ...
100  * </spectrum>
101  * </source>
102  *
103  * Optionally, a temporal model may be provided using the following
104  * syntax
105  *
106  * <source name="..." type="..." instrument="...">
107  * <spectrum type="...">
108  * ...
109  * </spectrum>
110  * <temporalModel type="...">
111  * ...
112  * </temporalModel>
113  * </source>
114  *
115  * If no temporal component is found a constant model is assumed.
116  ***************************************************************************/
118  GModelData(xml)
119 {
120  // Initialise members
121  init_members();
122 
123  // Read XML
124  read(xml);
125 
126  // Set parameter pointers
127  set_pointers();
128 
129  // Return
130  return;
131 }
132 
133 
134 /***********************************************************************//**
135  * @brief Construct from spectral component
136  *
137  * @param[in] spectral Spectral model component.
138  *
139  * Constructs a CTA instrumental background model from a spectral
140  * model component. The temporal component is assumed to be constant.
141  * Please refer to the classe GModelSpectral to learn more about the
142  * definition of the spectral components.
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 Copy constructor
167  *
168  * @param[in] bgd CTA instrument background model.
169  ***************************************************************************/
171  GModelData(bgd)
172 {
173  // Initialise class members
174  init_members();
175 
176  // Copy members
177  copy_members(bgd);
178 
179  // Return
180  return;
181 }
182 
183 
184 /***********************************************************************//**
185  * @brief Destructor
186  ***************************************************************************/
188 {
189  // Free members
190  free_members();
191 
192  // Return
193  return;
194 }
195 
196 
197 /*==========================================================================
198  = =
199  = Operators =
200  = =
201  ==========================================================================*/
202 
203 /***********************************************************************//**
204  * @brief Assignment operator
205  *
206  * @param[in] bgd CTA instrument background model.
207  * @return CTA instrument background model.
208  ***************************************************************************/
210 {
211  // Execute only if object is not identical
212  if (this != &bgd) {
213 
214  // Copy base class members
215  this->GModelData::operator=(bgd);
216 
217  // Free members
218  free_members();
219 
220  // Initialise private members
221  init_members();
222 
223  // Copy members
224  copy_members(bgd);
225 
226  } // endif: object was not identical
227 
228  // Return this object
229  return *this;
230 }
231 
232 
233 /*==========================================================================
234  = =
235  = Public methods =
236  = =
237  ==========================================================================*/
238 
239 /***********************************************************************//**
240  * @brief Clear CTA instrument background model
241  *
242  * This method properly resets the CTA instrument background model to an
243  * initial state.
244  ***************************************************************************/
246 {
247  // Free class members (base and derived classes, derived class first)
248  free_members();
249  this->GModelData::free_members();
250 
251  // Initialise members
252  this->GModelData::init_members();
253  init_members();
254 
255  // Return
256  return;
257 }
258 
259 
260 /***********************************************************************//**
261  * @brief Clone CTA instrument background model
262  *
263  * @return Pointer to deep copy of CTA instrument background model.
264  ***************************************************************************/
266 {
267  return new GCTAModelIrfBackground(*this);
268 }
269 
270 
271 /***********************************************************************//**
272  * @brief Set spectral model component
273  *
274  * @param[in] spectral Pointer to spectral model component.
275  *
276  * Sets the spectral model component of the model.
277  ***************************************************************************/
279 {
280  // Free spectral model component only if it differs from current
281  // component. This prevents unintential deallocation of the argument
282  if ((m_spectral != NULL) && (m_spectral != spectral)) {
283  delete m_spectral;
284  }
285 
286  // Clone spectral model component if it exists, otherwise set pointer
287  // to NULL
288  m_spectral = (spectral != NULL) ? spectral->clone() : NULL;
289 
290  // Set parameter pointers
291  set_pointers();
292 
293  // Return
294  return;
295 }
296 
297 
298 /***********************************************************************//**
299  * @brief Set temporal model component
300  *
301  * @param[in] temporal Pointer to temporal model component.
302  *
303  * Sets the temporal model component of the model.
304  ***************************************************************************/
306 {
307  // Free temporal model component only if it differs from current
308  // component. This prevents unintential deallocation of the argument
309  if ((m_temporal != NULL) && (m_temporal != temporal)) {
310  delete m_temporal;
311  }
312 
313  // Clone temporal model component if it exists, otherwise set pointer
314  // to NULL
315  m_temporal = (temporal != NULL) ? temporal->clone() : NULL;
316 
317  // Set parameter pointers
318  set_pointers();
319 
320  // Return
321  return;
322 }
323 
324 
325 /***********************************************************************//**
326  * @brief Return background rate in units of events MeV\f$^{-1}\f$
327  * s\f$^{-1}\f$ sr\f$^{-1}\f$
328  *
329  * @param[in] event Observed event.
330  * @param[in] obs Observation.
331  * @param[in] gradients Compute gradients?
332  * @return Background rate (events MeV\f$^{-1}\f$ s\f$^{-1}\f$ sr\f$^{-1}\f$).
333  *
334  * The method returns a real rate, defined by the number of counts per MeV,
335  * steradian and ontime.
336  *
337  * If the @p gradients flag is true the method will also set the parameter
338  * gradients of the model parameters.
339  ***************************************************************************/
341  const GObservation& obs,
342  const bool& gradients) const
343 {
344  // Get reference on CTA pointing and effective area response from
345  // observation and reference on CTA instrument direction from event
346  const GCTAPointing& pnt = gammalib::cta_pnt(G_EVAL, obs);
347  const GCTABackground& bgd = gammalib::cta_rsp_bkg(G_EVAL, obs);
348  const GCTAInstDir& dir = gammalib::cta_dir(G_EVAL, event);
349 
350  // Set DETX and DETY in instrument direction
351  GCTAInstDir inst_dir = pnt.instdir(dir.dir());
352 
353  // Evaluate function
354  double logE = event.energy().log10TeV();
355  double spat = bgd(logE, inst_dir.detx(), inst_dir.dety());
356  double spec = (spectral() != NULL)
357  ? spectral()->eval(event.energy(), event.time(), gradients)
358  : 1.0;
359  double temp = (temporal() != NULL)
360  ? temporal()->eval(event.time(), gradients) : 1.0;
361 
362  // Compute value
363  double value = spat * spec * temp;
364 
365  // If gradients were requested then multiply factors to spectral and
366  // temporal gradients
367  if (gradients) {
368 
369  // Multiply factors to spectral gradients
370  if (spectral() != NULL) {
371  double fact = spat * temp;
372  if (fact != 1.0) {
373  for (int i = 0; i < spectral()->size(); ++i)
374  (*spectral())[i].factor_gradient((*spectral())[i].factor_gradient() * fact);
375  }
376  }
377 
378  // Multiply factors to temporal gradients
379  if (temporal() != NULL) {
380  double fact = spat * spec;
381  if (fact != 1.0) {
382  for (int i = 0; i < temporal()->size(); ++i)
383  (*temporal())[i].factor_gradient((*temporal())[i].factor_gradient() * fact);
384  }
385  }
386 
387  } // endif: gradients were requested
388 
389  // Return value
390  return value;
391 }
392 
393 
394 /***********************************************************************//**
395  * @brief Return spatially integrated background rate in units of
396  * events MeV\f$^{-1}\f$ s\f$^{-1}\f$
397  *
398  * @param[in] obsEng Measured event energy.
399  * @param[in] obsTime Measured event time.
400  * @param[in] obs Observation.
401  * @return Spatially integrated background rate
402  * (events MeV\f$^{-1}\f$ s\f$^{-1}\f$)
403  *
404  * Spatially integrates the instrumental background model for a given
405  * measured event energy and event time. The method returns a real rate,
406  * defined as the number of counts per MeV and ontime.
407  ***************************************************************************/
409  const GTime& obsTime,
410  const GObservation& obs) const
411 {
412  // Set number of iterations for Romberg integration.
413  static const int iter_theta = 6;
414  static const int iter_phi = 6;
415 
416  // Initialise result
417  double npred = 0.0;
418  bool has_npred = false;
419 
420  // Build unique identifier
421  std::string id = obs.instrument() + "::" + obs.id();
422 
423  // Check if Npred value is already in cache
424  #if defined(G_USE_NPRED_CACHE)
425  if (!m_npred_names.empty()) {
426 
427  // Search for unique identifier, and if found, recover Npred value
428  // and break
429  for (int i = 0; i < m_npred_names.size(); ++i) {
430  if (m_npred_names[i] == id && m_npred_energies[i] == obsEng) {
431  npred = m_npred_values[i];
432  has_npred = true;
433  #if defined(G_DEBUG_NPRED)
434  std::cout << "GCTAModelIrfBackground::npred:";
435  std::cout << " cache=" << i;
436  std::cout << " npred=" << npred << std::endl;
437  #endif
438  break;
439  }
440  }
441 
442  } // endif: there were values in the Npred cache
443  #endif
444 
445  // Continue only if no Npred cache value has been found
446  if (!has_npred) {
447 
448  // Evaluate only if model is valid
449  if (valid_model()) {
450 
451  // Get reference on CTA pointing, background response and event
452  // list from observation
453  const GCTAPointing& pnt = gammalib::cta_pnt(G_NPRED, obs);
454  const GCTABackground& bgd = gammalib::cta_rsp_bkg(G_NPRED, obs);
455  const GCTAEventList& events = gammalib::cta_event_list(G_NPRED, obs);
456 
457  // Get instrument direction of RoI centre
458  GCTAInstDir roi_centre = pnt.instdir(events.roi().centre().dir());
459 
460  // Get ROI radius in radians
461  double roi_radius = events.roi().radius() * gammalib::deg2rad;
462 
463  // Get log10 of energy in TeV
464  double logE = obsEng.log10TeV();
465 
466  // Setup integration function
468  logE,
469  roi_centre,
470  iter_phi);
471 
472  // Setup integration
473  GIntegral integral(&integrand);
474 
475  // Set fixed number of iterations
476  integral.fixed_iter(iter_theta);
477 
478  // Spatially integrate radial component
479  npred = integral.romberg(0.0, roi_radius);
480 
481  // Store result in Npred cache
482  #if defined(G_USE_NPRED_CACHE)
483  m_npred_names.push_back(id);
484  m_npred_energies.push_back(obsEng);
485  m_npred_times.push_back(obsTime);
486  m_npred_values.push_back(npred);
487  #endif
488 
489  // Debug: Check for NaN
490  #if defined(G_NAN_CHECK)
491  if (gammalib::is_notanumber(npred) || gammalib::is_infinite(npred)) {
492  std::string origin = "GCTAModelIrfBackground::npred";
493  std::string message = " NaN/Inf encountered (npred=" +
494  gammalib::str(npred) + ", roi_radius=" +
495  gammalib::str(roi_radius) + ")";
496  gammalib::warning(origin, message);
497  }
498  #endif
499 
500  } // endif: model was valid
501 
502  } // endif: Npred computation required
503 
504  // Multiply in spectral and temporal components
505  npred *= spectral()->eval(obsEng, obsTime);
506  npred *= temporal()->eval(obsTime);
507 
508  // Return Npred
509  return npred;
510 }
511 
512 
513 /***********************************************************************//**
514  * @brief Return simulated list of events
515  *
516  * @param[in] obs Observation.
517  * @param[in] ran Random number generator.
518  * @return Pointer to list of simulated events (needs to be de-allocated by
519  * client)
520  *
521  * Draws a sample of events from the background model using a Monte
522  * Carlo simulation. The region of interest, the energy boundaries and the
523  * good time interval for the sampling will be extracted from the observation
524  * argument that is passed to the method. The method also requires a random
525  * number generator of type GRan which is passed by reference, hence the
526  * state of the random number generator will be changed by the method.
527  *
528  * For each event in the returned event list, the sky direction, the nominal
529  * coordinates (DETX and DETY), the energy and the time will be set.
530  ***************************************************************************/
532 {
533  // Initialise new event list
534  GCTAEventList* list = new GCTAEventList;
535 
536  // Continue only if model is valid)
537  if (valid_model()) {
538 
539  // Get reference on CTA pointing, background response and event list
540  // from observation
541  const GCTAPointing& pnt = gammalib::cta_pnt(G_MC, obs);
542  const GCTABackground& bgd = gammalib::cta_rsp_bkg(G_MC, obs);
543  const GCTAEventList& events = gammalib::cta_event_list(G_MC, obs);
544 
545  // Get simulation region
546  const GCTARoi& roi = events.roi();
547  const GEbounds& ebounds = events.ebounds();
548  const GGti& gti = events.gti();
549 
550  // Set simulation region for result event list
551  list->roi(roi);
552  list->ebounds(ebounds);
553  list->gti(gti);
554 
555  // Create a spectral model that combines the information from the
556  // background information and the spectrum provided by the model
558  for (int i = 0; i < spectral.nodes(); ++i) {
559  GEnergy energy = spectral.energy(i);
560  double intensity = spectral.intensity(i);
561  double norm = m_spectral->eval(energy);
562  spectral.intensity(i, norm*intensity);
563  }
564 
565  // Loop over all energy boundaries
566  for (int ieng = 0; ieng < ebounds.size(); ++ieng) {
567 
568  // Compute the background rate in model within the energy
569  // boundaries from spectral component (units: cts/s).
570  // Note that the time here is ontime. Deadtime correction will
571  // be done later.
572  double rate = spectral.flux(ebounds.emin(ieng), ebounds.emax(ieng));
573 
574  // Debug option: dump rate
575  #if defined(G_DUMP_MC)
576  std::cout << "GCTAModelIrfBackground::mc(\"" << name() << "\": ";
577  std::cout << "rate=" << rate << " cts/s)" << std::endl;
578  #endif
579 
580  // If the rate is not positive then skip this energy bins
581  if (rate <= 0.0) {
582  continue;
583  }
584 
585  // Loop over all good time intervals
586  for (int itime = 0; itime < gti.size(); ++itime) {
587 
588  // Get Monte Carlo event arrival times from temporal model
589  GTimes times = m_temporal->mc(rate,
590  gti.tstart(itime),
591  gti.tstop(itime),
592  ran);
593 
594  // Get number of events
595  int n_events = times.size();
596 
597  // Reserve space for events
598  if (n_events > 0) {
599  list->reserve(n_events);
600  }
601 
602  // Debug option: provide number of times and initialize
603  // statisics
604  #if defined(G_DUMP_MC)
605  std::cout << " Interval " << itime;
606  std::cout << " events=" << n_events << std::endl;
607  int n_killed_by_roi = 0;
608  #endif
609 
610  // Loop over events
611  for (int i = 0; i < n_events; ++i) {
612 
613  // Get Monte Carlo event energy from spectral model
614  GEnergy energy = spectral.mc(ebounds.emin(ieng),
615  ebounds.emax(ieng),
616  times[i],
617  ran);
618 
619  // Get Monte Carlo event direction from spatial model.
620  // This only will set the DETX and DETY coordinates.
621  GCTAInstDir instdir = bgd.mc(energy, times[i], ran);
622 
623  // Derive sky direction from instrument coordinates
624  GSkyDir skydir = pnt.skydir(instdir);
625 
626  // Set sky direction in GCTAInstDir object
627  instdir.dir(skydir);
628 
629  // Allocate event
630  GCTAEventAtom event;
631 
632  // Set event attributes
633  event.dir(instdir);
634  event.energy(energy);
635  event.time(times[i]);
636 
637  // Append event to list if it falls in ROI
638  if (events.roi().contains(event)) {
639  list->append(event);
640  }
641  #if defined(G_DUMP_MC)
642  else {
643  n_killed_by_roi++;
644  }
645  #endif
646 
647  } // endfor: looped over all events
648 
649  // Debug option: provide statisics
650  #if defined(G_DUMP_MC)
651  std::cout << " Killed by ROI=";
652  std::cout << n_killed_by_roi << std::endl;
653  #endif
654 
655  } // endfor: looped over all GTIs
656 
657  } // endfor: looped over all energy boundaries
658 
659  } // endif: model was valid
660 
661  // Return
662  return list;
663 }
664 
665 
666 /***********************************************************************//**
667  * @brief Read CTA instrument background model from XML element
668  *
669  * @param[in] xml XML element.
670  *
671  * Set up CTA instrument background model from the information provided by
672  * an XML element. The XML element is expected to have the following
673  * structure
674  *
675  * <source name="..." type="..." instrument="...">
676  * <spectrum type="...">
677  * ...
678  * </spectrum>
679  * </source>
680  *
681  * Optionally, a temporal model may be provided using the following
682  * syntax
683  *
684  * <source name="..." type="..." instrument="...">
685  * <spectrum type="...">
686  * ...
687  * </spectrum>
688  * <temporalModel type="...">
689  * ...
690  * </temporalModel>
691  * </source>
692  *
693  * If no temporal component is found a constant model is assumed.
694  ***************************************************************************/
696 {
697  // Clear model
698  clear();
699 
700  // Initialise XML elements
701  const GXmlElement* spectral = NULL;
702  const GXmlElement* temporal = NULL;
703 
704  // Get pointer on spectrum
705  spectral = xml.element("spectrum", 0);
706 
707  // Extract spectral model
708  m_spectral = xml_spectral(*spectral);
709 
710  // Optionally get temporal model
711  if (xml.elements("temporalModel")) {
712  temporal = xml.element("temporalModel", 0);
713  m_temporal = xml_temporal(*temporal);
714  }
715  else {
716  GModelTemporalConst constant;
717  m_temporal = constant.clone();
718  }
719 
720  // Read model attributes
721  read_attributes(xml);
722 
723  // Set parameter pointers
724  set_pointers();
725 
726  // Return
727  return;
728 }
729 
730 
731 /***********************************************************************//**
732  * @brief Write CTA instrument background model into XML element
733  *
734  * @param[in] xml XML element.
735  *
736  * Write CTA instrument background model information into an XML element.
737  * The XML element will have the following structure
738  *
739  * <source name="..." type="..." instrument="...">
740  * <spectrum type="...">
741  * ...
742  * </spectrum>
743  * </source>
744  *
745  * If the model contains a non-constant temporal model, the temporal
746  * component will also be written following the syntax
747  *
748  * <source name="..." type="..." instrument="...">
749  * <spectrum type="...">
750  * ...
751  * </spectrum>
752  * <temporalModel type="...">
753  * ...
754  * </temporalModel>
755  * </source>
756  *
757  * If no temporal component is found a constant model is assumed.
758  ***************************************************************************/
760 {
761  // Initialise pointer on source
762  GXmlElement* src = NULL;
763 
764  // Search corresponding source
765  int n = xml.elements("source");
766  for (int k = 0; k < n; ++k) {
767  GXmlElement* element = xml.element("source", k);
768  if (element->attribute("name") == name()) {
769  src = element;
770  break;
771  }
772  }
773 
774  // If we have a temporal model that is either not a constant, or a
775  // constant with a normalization value that differs from 1.0 then
776  // write the temporal component into the XML element. This logic
777  // assures compatibility with the Fermi/LAT format as this format
778  // does not handle temporal components.
779  bool write_temporal = ((m_temporal != NULL) &&
780  (m_temporal->type() != "Constant" ||
781  (*m_temporal)[0].value() != 1.0));
782 
783  // If no source with corresponding name was found then append one
784  if (src == NULL) {
785  src = xml.append("source");
786  if (spectral() != NULL) src->append(GXmlElement("spectrum"));
787  if (write_temporal) src->append(GXmlElement("temporalModel"));
788  }
789 
790  // Write spectral model
791  if (spectral() != NULL) {
792  GXmlElement* spec = src->element("spectrum", 0);
793  spectral()->write(*spec);
794  }
795 
796  // Optionally write temporal model
797  if (write_temporal) {
798  if (dynamic_cast<GModelTemporalConst*>(temporal()) == NULL) {
799  GXmlElement* temp = src->element("temporalModel", 0);
800  temporal()->write(*temp);
801  }
802  }
803 
804  // Write model attributes
805  write_attributes(*src);
806 
807  // Return
808  return;
809 }
810 
811 
812 /***********************************************************************//**
813  * @brief Print CTA instrument background model information
814  *
815  * @param[in] chatter Chattiness.
816  * @return String containing CTA instrument background model information.
817  ***************************************************************************/
818 std::string GCTAModelIrfBackground::print(const GChatter& chatter) const
819 {
820  // Initialise result string
821  std::string result;
822 
823  // Continue only if chatter is not silent
824  if (chatter != SILENT) {
825 
826  // Append header
827  result.append("=== GCTAModelIrfBackground ===");
828 
829  // Determine number of parameters per type
830  int n_spectral = (spectral() != NULL) ? spectral()->size() : 0;
831  int n_temporal = (temporal() != NULL) ? temporal()->size() : 0;
832 
833  // Append attributes
834  result.append("\n"+print_attributes());
835 
836  // Append model type
837  result.append("\n"+gammalib::parformat("Model type"));
838  if (n_spectral > 0) {
839  result.append("\""+spectral()->type()+"\"");
840  if (n_temporal > 0) {
841  result.append(" * ");
842  }
843  }
844  if (n_temporal > 0) {
845  result.append("\""+temporal()->type()+"\"");
846  }
847 
848  // Append parameters
849  result.append("\n"+gammalib::parformat("Number of parameters") +
850  gammalib::str(size()));
851  result.append("\n"+gammalib::parformat("Number of spectral par's") +
852  gammalib::str(n_spectral));
853  for (int i = 0; i < n_spectral; ++i) {
854  result.append("\n"+(*spectral())[i].print());
855  }
856  result.append("\n"+gammalib::parformat("Number of temporal par's") +
857  gammalib::str(n_temporal));
858  for (int i = 0; i < n_temporal; ++i) {
859  result.append("\n"+(*temporal())[i].print());
860  }
861 
862  } // endif: chatter was not silent
863 
864  // Return result
865  return result;
866 }
867 
868 
869 /*==========================================================================
870  = =
871  = Private methods =
872  = =
873  ==========================================================================*/
874 
875 /***********************************************************************//**
876  * @brief Initialise class members
877  ***************************************************************************/
879 {
880  // Initialise members
881  m_spectral = NULL;
882  m_temporal = NULL;
883 
884  // Initialise Npred cache
885  m_npred_names.clear();
886  m_npred_energies.clear();
887  m_npred_times.clear();
888  m_npred_values.clear();
889 
890  // Return
891  return;
892 }
893 
894 
895 /***********************************************************************//**
896  * @brief Copy class members
897  *
898  * @param[in] bgd CTA background model.
899  ***************************************************************************/
901 {
902  // Copy cache
907 
908  // Clone spectral and temporal model components
909  m_spectral = (bgd.m_spectral != NULL) ? bgd.m_spectral->clone() : NULL;
910  m_temporal = (bgd.m_temporal != NULL) ? bgd.m_temporal->clone() : NULL;
911 
912  // Set parameter pointers
913  set_pointers();
914 
915  // Return
916  return;
917 }
918 
919 
920 /***********************************************************************//**
921  * @brief Delete class members
922  ***************************************************************************/
924 {
925  // Free memory
926  if (m_spectral != NULL) delete m_spectral;
927  if (m_temporal != NULL) delete m_temporal;
928 
929  // Signal free pointers
930  m_spectral = NULL;
931  m_temporal = NULL;
932 
933  // Return
934  return;
935 }
936 
937 
938 /***********************************************************************//**
939  * @brief Set pointers
940  *
941  * Set pointers to all model parameters. The pointers are stored in a vector
942  * that is member of the GModelData base class.
943  ***************************************************************************/
945 {
946  // Clear parameters
947  m_pars.clear();
948 
949  // Determine the number of parameters
950  int n_spectral = (spectral() != NULL) ? spectral()->size() : 0;
951  int n_temporal = (temporal() != NULL) ? temporal()->size() : 0;
952  int n_pars = n_spectral + n_temporal;
953 
954  // Continue only if there are parameters
955  if (n_pars > 0) {
956 
957  // Gather spectral parameters
958  for (int i = 0; i < n_spectral; ++i) {
959  m_pars.push_back(&((*spectral())[i]));
960  }
961 
962  // Gather temporal parameters
963  for (int i = 0; i < n_temporal; ++i) {
964  m_pars.push_back(&((*temporal())[i]));
965  }
966 
967  }
968 
969  // Return
970  return;
971 }
972 
973 
974 /***********************************************************************//**
975  * @brief Verifies if model has all components
976  *
977  * Returns 'true' if models has a spectral and a temporal component.
978  * Otherwise returns 'false'.
979  ***************************************************************************/
981 {
982  // Set result
983  bool result = ((spectral() != NULL) && (temporal() != NULL));
984 
985  // Return result
986  return result;
987 }
988 
989 
990 /***********************************************************************//**
991  * @brief Return pointer to spectral model from XML element
992  *
993  * @param[in] spectral XML element.
994  * @return Pointer to spectral model.
995  *
996  * Returns pointer to spectral model that is defined in an XML element.
997  ***************************************************************************/
999 {
1000  // Get spectral model
1001  GModelSpectralRegistry registry;
1002  GModelSpectral* ptr = registry.alloc(spectral);
1003 
1004  // Return pointer
1005  return ptr;
1006 }
1007 
1008 
1009 /***********************************************************************//**
1010  * @brief Return pointer to temporal model from XML element
1011  *
1012  * @param[in] temporal XML element.
1013  * @return Pointer to temporal model.
1014  *
1015  * Returns pointer to temporal model that is defined in an XML element.
1016  ***************************************************************************/
1018 {
1019  // Get temporal model
1020  GModelTemporalRegistry registry;
1021  GModelTemporal* ptr = registry.alloc(temporal);
1022 
1023  // Return pointer
1024  return ptr;
1025 }
1026 
1027 
1028 /***********************************************************************//**
1029  * @brief Kernel for offset angle integration of background model
1030  *
1031  * @param[in] theta Offset angle from ROI centre (radians).
1032  *
1033  * Computes
1034  *
1035  * \f[
1036  * K(\rho | E, t) = \sin \theta \times
1037  * \int_{0}^{2\pi}
1038  * B(\theta,\phi | E, t) d\phi
1039  * \f]
1040  *
1041  * where \f$B(\theta,\phi | E, t)\f$ is the background model for a specific
1042  * observed energy \f$E\f$ and time \f$t\f$.
1043  ***************************************************************************/
1045 {
1046  // Initialise value
1047  double value = 0.0;
1048 
1049  // Continue only if offset angle is positive
1050  if (theta > 0.0) {
1051 
1052  // Setup phi integration kernel
1054  m_logE,
1055  m_roi_centre,
1056  theta);
1057 
1058  // Setup integration
1059  GIntegral integral(&integrand);
1060 
1061  // Set fixed number of iterations
1062  integral.fixed_iter(m_iter);
1063 
1064  // Integrate over phi
1065  value = integral.romberg(0.0, gammalib::twopi) * std::sin(theta);
1066 
1067  // Debug: Check for NaN
1068  #if defined(G_NAN_CHECK)
1069  if (gammalib::is_notanumber(value) || gammalib::is_infinite(value)) {
1070  std::string origin = "GCTAModelIrfBackground::npred_roi_kern_theta::eval"
1071  "(" + gammalib::str(theta) + ")";
1072  std::string message = " NaN/Inf encountered (value=" +
1073  gammalib::str(value) + ")";
1074  gammalib::warning(origin, message);
1075  }
1076  #endif
1077 
1078  } // endif: offset angle was positive
1079 
1080  // Return value
1081  return value;
1082 }
1083 
1084 
1085 /***********************************************************************//**
1086  * @brief Kernel for azimuth angle integration of background model
1087  *
1088  * @param[in] phi Azimuth angle around ROI centre (radians).
1089  *
1090  * Computes
1091  *
1092  * \f[
1093  * B(\theta, \phi | E, t)
1094  * \f]
1095  *
1096  * using
1097  *
1098  * \f[ {\rm detx} = \theta \cos \phi \f]
1099  * \f[ {\rm dety} = \theta \sin \phi \f]
1100  *
1101  * @todo Verify correct orientation of detx and dety with respect to phi
1102  ***************************************************************************/
1104 {
1105  // Compute detx and dety in radians
1106  double detx = m_roi_centre.detx();
1107  double dety = m_roi_centre.dety();
1108  if (m_theta > 0.0 ) {
1109  detx += m_theta * std::cos(phi);
1110  dety += m_theta * std::sin(phi);
1111  }
1112 
1113  // Get background value
1114  double value = (*m_bgd)(m_logE, detx, dety);
1115 
1116  // Debug: Check for NaN
1117  #if defined(G_NAN_CHECK)
1118  if (gammalib::is_notanumber(value) || gammalib::is_infinite(value)) {
1119  std::string origin = "GCTAModelIrfBackground::npred_roi_kern_phi::eval"
1120  "(" + gammalib::str(phi) + ")";
1121  std::string message = " NaN/Inf encountered (value=" +
1122  gammalib::str(value) + ", detx=" +
1123  gammalib::str(detx) + ", dety=" +
1124  gammalib::str(dety) + ")";
1125  gammalib::warning(origin, message);
1126  }
1127  #endif
1128 
1129  // Return Npred
1130  return value;
1131 }
void copy_members(const GCTAModelIrfBackground &bgd)
Copy class members.
virtual void write(GXmlElement &xml) const =0
CTA background model base class definition.
const GCTABackground * m_bgd
Pointer to background.
Constant temporal model class interface definition.
std::vector< GEnergy > m_npred_energies
Model energy.
virtual double eval(const GEvent &event, const GObservation &obs, const bool &gradients=false) const
Return background rate in units of events MeV s sr .
GModelTemporal * m_temporal
Temporal model.
std::vector< GTime > m_npred_times
Model time.
virtual GCTAEventList * mc(const GObservation &obs, GRan &ran) const
Return simulated list of events.
#define G_EVAL
virtual GTimes mc(const double &rate, const GTime &tmin, const GTime &tmax, GRan &ran) const =0
virtual GEnergy mc(const GEnergy &emin, const GEnergy &emax, const GTime &time, GRan &ran) const =0
double norm(const GVector &vector)
Computes vector norm.
Definition: GVector.cpp:821
double romberg(std::vector< double > bounds, const int &order=5)
Perform Romberg integration.
Definition: GIntegral.cpp:380
void detx(const double &x)
Set DETX coordinate (in radians)
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
CTA IRF background model class.
CTA IRF background model class definition.
Abstract temporal model base class.
virtual void roi(const GRoi &roi)
Set ROI.
GSkyDir skydir(const GCTAInstDir &instdir) const
Get sky direction direction from instrument direction.
virtual std::string type(void) const
Return data model type.
int size(void) const
Return number of parameters.
void ebounds(const GEbounds &ebounds)
Set energy boundaries.
Definition: GEvents.cpp:136
virtual GCTAModelIrfBackground * clone(void) const
Clone CTA instrument background model.
int size(void) const
Return number of energy boundaries.
Definition: GEbounds.hpp:157
GVector cos(const GVector &vector)
Computes cosine of vector elements.
Definition: GVector.cpp:1100
virtual double eval(const GEnergy &srcEng, const GTime &srcTime=GTime(), const bool &gradients=false) const =0
void append(const GCTAEventAtom &event)
Append event to event list.
GCTAModelIrfBackground(void)
Void constructor.
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.
XML element node class.
Definition: GXmlElement.hpp:47
void gti(const GGti &gti)
Set Good Time Intervals.
Definition: GEvents.cpp:154
void free_members(void)
Delete class members.
Definition: GModelData.cpp:181
Spectral model registry class definition.
GModelTemporal * xml_temporal(const GXmlElement &temporal) const
Return pointer to temporal model from XML element.
virtual void write(GXmlElement &xml) const =0
double log10TeV(void) const
Return log10 of energy in TeV.
Definition: GEnergy.cpp:457
Random number generator class.
Definition: GRan.hpp:44
virtual const GModelSpectralNodes & spectrum(void) const =0
Spectral nodes model class definition.
bool valid_model(void) const
Verifies if model has all components.
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
double eval(const double &theta)
Kernel for offset angle integration of background model.
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.
void free_members(void)
Delete class members.
virtual const GTime & time(void) const =0
void init_members(void)
Initialise class members.
const GCTABackground & cta_rsp_bkg(const std::string &origin, const GObservation &obs)
Retrieve CTA background response from generic observation.
void dir(const GSkyDir &dir)
Set sky direction.
void set_pointers(void)
Set pointers.
int size(void) const
Return number of Good Time Intervals.
Definition: GGti.hpp:153
#define G_NPRED
Spectral nodes model class.
virtual std::string print(const GChatter &chatter=NORMAL) const
Print CTA instrument background model information.
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.
const double deg2rad
Definition: GMath.hpp:43
GCTAModelIrfBackground & operator=(const GCTAModelIrfBackground &bgd)
Assignment operator.
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.
const GEnergy & emin(void) const
Return minimum energy of all intervals.
Definition: GEbounds.hpp:181
virtual void write(GXmlElement &xml) const
Write CTA instrument background model into XML element.
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
virtual ~GCTAModelIrfBackground(void)
Destructor.
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 GModelData & operator=(const GModelData &model)
Assignment operator.
Definition: GModelData.cpp:120
GChatter
Definition: GTypemaps.hpp:33
std::vector< double > m_npred_values
Model values.
GModelSpectral * spectral(void) const
Return spectral model component.
Good Time Interval class.
Definition: GGti.hpp:62
virtual double eval(const GTime &srcTime, const bool &gradients=false) const =0
Abstract observation base class.
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
#define G_MC
CTA observation class interface definition.
virtual void clear(void)
Clear CTA instrument background model.
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.
GModelSpectral * xml_spectral(const GXmlElement &spectral) const
Return pointer to spectral model from XML element.
GModelTemporal * temporal(void) const
Return temporal model component.
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.
const GCTAModelCubeBackground g_cta_inst_background_seed
virtual void read(const GXmlElement &xml)
Read CTA instrument background model from XML element.
CTA event atom class.
GVector sin(const GVector &vector)
Computes sine of vector elements.
Definition: GVector.cpp:1226
double eval(const double &phi)
Kernel for azimuth angle integration of background model.
Abstract base class for the CTA background model.
void fixed_iter(const int &iter)
Set fixed number of iterations.
Definition: GIntegral.hpp:181
GModelSpectral * m_spectral
Spectral model.
virtual GModelTemporal * clone(void) const =0
Clones object.
void read_attributes(const GXmlElement &xml)
Read model attributes.
Definition: GModel.cpp:668
CTA instrument direction class.
Definition: GCTAInstDir.hpp:59
virtual GCTAInstDir mc(const GEnergy &energy, const GTime &time, GRan &ran) const =0
void dety(const double &y)
Set DETY coordinate (in radians)
const double twopi
Definition: GMath.hpp:36
const GEnergy & emax(void) const
Return maximum energy of all intervals.
Definition: GEbounds.hpp:193
virtual double flux(const GEnergy &emin, const GEnergy &emax) const =0
Model registry class definition.
std::vector< std::string > m_npred_names
Model names.
virtual GXmlNode * append(const GXmlNode &node)
Append XML child node.
Definition: GXmlNode.cpp:284
std::string parformat(const std::string &s, const int &indent=0)
Convert string in parameter format.
Definition: GTools.cpp:1022
virtual double npred(const GEnergy &obsEng, const GTime &obsTime, const GObservation &obs) const
Return spatially integrated background rate in units of events MeV s .
Integration class interface definition.
Sky direction class.
Definition: GSkyDir.hpp:62
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