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