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