GammaLib  1.7.0.dev
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
GCTASupport.cpp
Go to the documentation of this file.
1 /***************************************************************************
2  * GCTASupport.cpp - CTA support functions *
3  * ----------------------------------------------------------------------- *
4  * copyright (C) 2011-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 GCTASupport.cpp
23  * @brief Implementation of support function used by CTA classes
24  * @author Juergen Knoedlseder
25  */
26 
27 /* __ Includes ___________________________________________________________ */
28 #ifdef HAVE_CONFIG_H
29 #include <config.h>
30 #endif
31 #include <typeinfo>
32 #include <cmath>
33 #include <iostream>
34 #include <string>
35 #include "GCTASupport.hpp"
36 #include "GTools.hpp"
37 #include "GMath.hpp"
38 #include "GFitsHDU.hpp"
39 #include "GEbounds.hpp"
40 #include "GPhases.hpp"
41 #include "GResponse.hpp"
42 #include "GModelData.hpp"
43 #include "GModelSpectral.hpp"
44 #include "GModelTemporal.hpp"
45 #include "GCTAObservation.hpp"
46 #include "GCTAPointing.hpp"
47 #include "GCTAResponse.hpp"
48 #include "GCTAResponseIrf.hpp"
49 #include "GCTAResponseCube.hpp"
50 #include "GCTAAeff.hpp"
51 #include "GCTABackground.hpp"
52 #include "GCTAEventList.hpp"
53 #include "GCTARoi.hpp"
54 #include "GCTAInstDir.hpp"
55 #include "GCTAModelBackground.hpp"
60 
61 /* __ Method name definitions ____________________________________________ */
62 #define G_READ_DS_ROI "gammalib::read_ds_roi(GFitsHDU&)"
63 #define G_READ_DS_EBOUNDS "gammalib::read_ds_ebounds(GFitsHDU&)"
64 #define G_READ_DS_PHASE "gammalib::read_ds_phase(GFitsHDU&)"
65 #define G_READ_DS_GTI "gammalib::read_ds_gti_extname(GFitsHDU&)"
66 #define G_CTA_MODEL_NROI "gammalib::cta_model_nroi(GObservation*, "\
67  "GModelData*, GEnergy&, GTime&)"
68 
69 /* __ Coding definitions _________________________________________________ */
70 
71 /* __ Debug definitions __________________________________________________ */
72 
73 
74 /***********************************************************************//**
75  * @brief Extract ROI from data sub-space keywords
76  *
77  * @param[in] hdu FITS HDU
78  *
79  * @exception GException::invalid_value
80  * Invalid ROI data sub-space encountered
81  *
82  * Reads the ROI data sub-space keywords by searching for a DSTYPx keyword
83  * named "POS(RA,DEC)". The data sub-space information is expected to be in
84  * the format "CIRCLE(267.0208,-24.78,4.5)", where the 3 arguments are Right
85  * Ascension, Declination and radius in units of degrees. No detailed syntax
86  * checking is performed.
87  *
88  * If no ROI information has been found, an GCTARoi object with initial
89  * values will be returned.
90  ***************************************************************************/
92 {
93  // Initialise ROI
94  GCTARoi roi;
95 
96  // Get number of data sub-space keywords (default to 0 if keyword is
97  // not found)
98  int ndskeys = (hdu.has_card("NDSKEYS")) ? hdu.integer("NDSKEYS") : 0;
99 
100  // Loop over all data selection keys
101  for (int i = 1; i <= ndskeys; ++i) {
102 
103  // Set data sub-space key strings
104  std::string type_key = "DSTYP"+gammalib::str(i);
105  //std::string unit_key = "DSUNI"+gammalib::str(i);
106  std::string value_key = "DSVAL"+gammalib::str(i);
107 
108  // Continue only if type_key is found and if this key is POS(RA,DEC)
109  if (hdu.has_card(type_key) && hdu.string(type_key) == "POS(RA,DEC)") {
110 
111  // ...
112  //std::string unit = gammalib::toupper(hdu.string(unit_key));
113  std::string value = hdu.string(value_key);
114  std::string value_proc = gammalib::strip_chars(value, "CIRCLE(");
115  value_proc = gammalib::strip_chars(value_proc, ")");
116  std::vector<std::string> args = gammalib::split(value_proc, ",");
117  if (args.size() == 3) {
118 
119  // Extract sky direction and radius
120  double ra = gammalib::todouble(args[0]);
121  double dec = gammalib::todouble(args[1]);
122  double rad = gammalib::todouble(args[2]);
123  GSkyDir skydir;
124  skydir.radec_deg(ra, dec);
125 
126  // If the header contains pointing information then set a
127  // full instrument direction, including DETX and DETY
128  // coordinates.
129  if (hdu.has_card("RA_PNT") && hdu.has_card("DEC_PNT")) {
130  double ra_pnt = hdu.real("RA_PNT");
131  double dec_pnt = hdu.real("DEC_PNT");
132  GSkyDir dir_pnt;
133  dir_pnt.radec_deg(ra_pnt, dec_pnt);
134  GCTAPointing pnt(dir_pnt);
135  GCTAInstDir dir = pnt.instdir(skydir);
136  roi.centre(dir);
137  }
138 
139  // ... otherwise just set the sky direction
140  else {
141  GCTAInstDir dir(skydir);
142  roi.centre(dir);
143  }
144 
145  // Set RoI radius
146  roi.radius(rad);
147  }
148  else {
149  std::string msg = "Invalid acceptance cone value \""+value+
150  "\" encountered in data sub-space "
151  "key \""+value_key+"\".";
153  }
154 
155  } // endif: POS(RA,DEC) type found
156 
157  } // endfor: looped over data sub-space keys
158 
159  // Return roi
160  return roi;
161 }
162 
163 
164 /***********************************************************************//**
165  * @brief Read energy boundary data sub-space keywords
166  *
167  * @param[in] hdu FITS HDU
168  *
169  * @exception GException::invalid_value
170  * Invalid energy data sub-space encountered
171  *
172  * Reads the energy boundary data sub-space keywords by searching for a
173  * DSTYPx keyword named "ENERGY". The data sub-space information is expected
174  * to be in the format "200:50000", where the 2 arguments are the minimum
175  * and maximum energy. The energy unit is given by the keyword DSUNIx, which
176  * supports keV, MeV, GeV and TeV (case independent). No detailed syntax
177  * checking is performed.
178  ***************************************************************************/
180 {
181  // Initialise energy boundaries
182  GEbounds ebounds;
183 
184  // Get number of data sub-space keywords (default to 0 if keyword is
185  // not found)
186  int ndskeys = (hdu.has_card("NDSKEYS")) ? hdu.integer("NDSKEYS") : 0;
187 
188  // Loop over all data sub-space keys
189  for (int i = 1; i <= ndskeys; ++i) {
190 
191  // Set data sub-space key strings
192  std::string type_key = "DSTYP"+gammalib::str(i);
193  std::string unit_key = "DSUNI"+gammalib::str(i);
194  std::string value_key = "DSVAL"+gammalib::str(i);
195 
196  // Continue only if type_key is found and if this key is ENERGY
197  if (hdu.has_card(type_key) && hdu.string(type_key) == "ENERGY") {
198 
199  // Extract energy boundaries
200  std::string value = hdu.string(value_key);
201  std::string unit = hdu.string(unit_key);
202  std::vector<std::string> values = gammalib::split(value, ":");
203  if (values.size() == 2) {
204  double emin = gammalib::todouble(values[0]);
205  double emax = gammalib::todouble(values[1]);
206  GEnergy e_min(emin, unit);
207  GEnergy e_max(emax, unit);
208  ebounds.append(e_min, e_max);
209  }
210  else {
211  std::string msg = "Invalid energy value \""+value+
212  "\" encountered in data selection key \""+
213  value_key+"\"";
215  }
216 
217  } // endif: ENERGY type_key found
218 
219  } // endfor: looped over data selection keys
220 
221  // Return
222  return ebounds;
223 }
224 
225 
226 /***********************************************************************//**
227  * @brief Read phase boundary data sub-space keywords
228  *
229  * @param[in] hdu FITS HDU
230  * @return Phase intervals
231  *
232  * @exception GException::invalid_value
233  * Invalid phase data sub-space encountered
234  *
235  * Reads the phase boundary data sub-space keywords by searching for a
236  * DSTYPx keyword named "PHASE". The data sub-space information is expected
237  * to be in the format "0.1:0.3,0.5:0.7", where each subset of numbers
238  * represents the minimum and maximum phase. No detailed syntax
239  * checking is performed.
240  ***************************************************************************/
242 {
243  // Initialise phase intervals
244  GPhases phases;
245 
246  // Get number of data sub-space keywords (default to 0 if keyword is
247  // not found)
248  int ndskeys = (hdu.has_card("NDSKEYS")) ? hdu.integer("NDSKEYS") : 0;
249 
250  // Loop over all data sub-space keys
251  for (int i = 1; i <= ndskeys; ++i) {
252 
253  // Set data sub-space key strings
254  std::string type_key = "DSTYP"+gammalib::str(i);
255  std::string unit_key = "DSUNI"+gammalib::str(i);
256  std::string value_key = "DSVAL"+gammalib::str(i);
257 
258  // Continue only if type_key is found and if this key is PHASE
259  if (hdu.has_card(type_key) && hdu.string(type_key) == "PHASE") {
260 
261  // Extract phase boundaries
262  std::string value = hdu.string(value_key);
263  std::vector<std::string> intervals = gammalib::split(value, ",");
264  for (int j = 0; j < intervals.size(); ++j) {
265  std::vector<std::string> values = gammalib::split(intervals[j], ":");
266  if (values.size() == 2) {
267  double phmin = gammalib::todouble(values[0]);
268  double phmax = gammalib::todouble(values[1]);
269  phases.append(phmin, phmax);
270  }
271  else {
272  std::string msg = "Invalid phase value \""+value+
273  "\" encountered in data selection key \""+
274  value_key+"\"";
276  }
277 
278  } //endfor: looped over phase intervals
279 
280  } // endif: PHASE type_key found
281 
282  } // endfor: looped over data selection keys
283 
284  // Return
285  return phases;
286 }
287 
288 
289 /***********************************************************************//**
290  * @brief Return Good Time Intervals extension name from data sub-space
291  * keywords
292  *
293  * @param[in] hdu FITS HDU
294  * @return Good Time Interval extension name
295  *
296  * @exception GException::invalid_value
297  * Invalid Good Time Intervals data sub-space encountered
298  *
299  * Returns the name of the FITS extension that contains the Good Time
300  * Intervals by screening the data sub-space keywords that are present in
301  * the FITS header. The method searches for a DSTYPx keyword named "TIME"
302  * and a corresponding DSVALx keyword named "TABLE", and the extension name
303  * is extracted from the corresponding DSREFx keyword. Note that by
304  * convention the extension name is preceeded by a colon, which is stripped
305  * by this method.
306  ***************************************************************************/
308 {
309  // Initialise extension name
310  std::string extname;
311 
312  // Get number of data sub-space keys (default to 0 if "NDSKEYS" keyword
313  // is not found)
314  int ndskeys = (hdu.has_card("NDSKEYS")) ? hdu.integer("NDSKEYS") : 0;
315 
316  // Loop over all data sub-space keys
317  for (int i = 1; i <= ndskeys; ++i) {
318 
319  // Set data sub-space key strings
320  std::string type_key = "DSTYP"+gammalib::str(i);
321  std::string unit_key = "DSUNI"+gammalib::str(i);
322  std::string value_key = "DSVAL"+gammalib::str(i);
323  std::string ref_key = "DSREF"+gammalib::str(i);
324 
325  // Skip if DSTYPi keyword does not exist, or if it is not "TIME"
326  if (!hdu.has_card(type_key) || hdu.string(type_key) != "TIME") {
327  continue;
328  }
329 
330  // If DSVALi keyword does not exist or if it's value is not
331  // "TABLE" then throw an exception
332  if (!hdu.has_card(value_key)) {
333  std::string msg = "Keyword \""+value_key+"\" missing for data "
334  "sub-space selection of type "+type_key+
335  "=\"TIME\". Please correct FITS header.";
337  }
338  if (hdu.string(value_key) != "TABLE") {
339  std::string msg = "Cannot interpret keyword \""+value_key+"\" "
340  "value \""+hdu.string(value_key)+"\". Only "
341  "the value \"TABLE\" is supported.";
343  }
344 
345  // If DSREFi keyword does not exist then throw an exception
346  if (!hdu.has_card(ref_key)) {
347  std::string msg = "Keyword \""+ref_key+"\" missing for data "
348  "sub-space selection of type "+type_key+
349  "=\"TIME\". Please correct FITS header.";
351  }
352 
353  // Get extension name (strip any leading and trailing colons)
354  extname = gammalib::strip_whitespace(gammalib::strip_chars(hdu.string(ref_key), ":"));
355 
356  } // endfor: looped over data sub-space keywords
357 
358  // Return
359  return extname;
360 }
361 
362 
363 /***********************************************************************//**
364  * @brief Return extension name for GADF response table of given HDU class 4
365  *
366  * @param[in] fits FITS file.
367  * @param[in] hduclas4 HDU class 4.
368  * @return Extension name.
369  *
370  * Returns the extension name for GADF response table of given HDU class 4.
371  * If the response table is not found, an empty extension name is returned.
372  ***************************************************************************/
373 std::string gammalib::gadf_hduclas4(const GFits& fits,
374  const std::string& hduclas4)
375 {
376  // Initialise extension name
377  std::string extname;
378 
379  // Loop over all FITS HDUs
380  for (int i = 0; i < fits.size(); ++i) {
381 
382  // Get FITS HDU
383  const GFitsHDU* hdu = fits[i];
384 
385  // Skip HDU if it has no "HDUCLASS", "HDUCLAS1" and "HDUCLAS4" keywords
386  if (!(hdu->has_card("HDUCLASS") &&
387  hdu->has_card("HDUCLAS1") &&
388  hdu->has_card("HDUCLAS4"))) {
389  continue;
390  }
391 
392  // Skip HDU if "HDUCLASS" is not "GADF", "HDUCLAS1" is not response
393  // and "HDUCLAS4" is not hduclas4
394  if (!((gammalib::strip_whitespace(hdu->string("HDUCLASS")) == "GADF") &&
395  (gammalib::strip_whitespace(hdu->string("HDUCLAS1")) == "RESPONSE") &&
396  (gammalib::strip_whitespace(hdu->string("HDUCLAS4")) == hduclas4))) {
397  continue;
398  }
399 
400  // Extract extension name and break
401  extname = hdu->extname();
402 
403  } // endfor: loop over headers
404 
405  // Return
406  return extname;
407 }
408 
409 
410 /***********************************************************************//**
411  * @brief Determine number of radial Romberg iterations
412  *
413  * @param[in] rho_max Maximum radial offset (radians).
414  * @param[in] resolution Requested angular resolution (radians).
415  * @param[in] iter_min Minimum number of iterations.
416  * @param[in] iter_max Maximum number of iterations.
417  * @return Number of radial Romberg iterations.
418  *
419  * Determines the number of radial Romberg iterations using the formula
420  *
421  * \f[
422  * iter = \log_2 \left( \frac{\rho_{\rm max}}{resolution} \right) + 1
423  * \f]
424  *
425  * where
426  * \f$\rho_{\rm max}\f$ is the maximum radial offset and
427  * \f$resolution\f$ is the required angular resolution.
428  *
429  * The result will be constrained to the interval [@p iter_min,@p iter_max].
430  ***************************************************************************/
431 int gammalib::iter_rho(const double& rho_max,
432  const double& resolution,
433  const int& iter_min,
434  const int& iter_max)
435 {
436  // Initialise number of iterations
437  int iter = iter_min;
438 
439  // Continue only if radial offset and resolution are positive
440  if ((rho_max > 0.0) && (resolution > 0.0)) {
441  double arg = rho_max / resolution;
442  iter = int(std::log(arg) * gammalib::inv_ln2) + 1;
443  if (iter < iter_min) {
444  iter = iter_min;
445  }
446  else if (iter > iter_max) {
447  iter = iter_max;
448  }
449  }
450 
451  // Return number of iterations
452  return iter;
453 }
454 
455 
456 /***********************************************************************//**
457  * @brief Determine number of azimuthal Romberg iterations
458  *
459  * @param[in] rho Radial offset (radians).
460  * @param[in] resolution Requested angular resolution (radians).
461  * @param[in] iter_min Minimum number of iterations.
462  * @param[in] iter_max Maximum number of iterations.
463  * @return Number of azimuthal Romberg iterations.
464  *
465  * Determines the number of azimuthal Romberg iterations using the formula
466  *
467  * \f[
468  * iter = \log_2 \left( \frac{2\pi \rho}{resolution} \right) + 1
469  * \f]
470  *
471  * where
472  * \f$\rho\f$ is the radial offset and
473  * \f$resolution\f$ is the required angular resolution.
474  *
475  * The result will be constrained to the interval [@p iter_min,@p iter_max].
476  ***************************************************************************/
477 int gammalib::iter_phi(const double& rho,
478  const double& resolution,
479  const int& iter_min,
480  const int& iter_max)
481 {
482  // Initialise number of iterations
483  int iter = iter_min;
484 
485  // Continue only if radial offset and resolution are positive
486  if ((rho > 0.0) && (resolution > 0.0)) {
487  double arg = gammalib::twopi * rho / resolution;
488  iter = int(std::log(arg) * gammalib::inv_ln2) + 1;
489  if (iter < iter_min) {
490  iter = iter_min;
491  }
492  else if (iter > iter_max) {
493  iter = iter_max;
494  }
495  }
496 
497  // Return number of iterations
498  return iter;
499 }
500 
501 
502 /***********************************************************************//**
503  * @brief Retrieve CTA observation from generic observation
504  *
505  * @param[in] origin Method asking for pointer retrieval.
506  * @param[in] obs Generic observation.
507  * @return Reference to CTA observation.
508  *
509  * @exception GException::invalid_argument
510  * Observation @p obs is not a CTA observations.
511  *
512  * Dynamically casts generic observation into a CTA observation. If the
513  * generic observation is not a CTA observation, an exception is thrown.
514  ***************************************************************************/
515 const GCTAObservation& gammalib::cta_obs(const std::string& origin,
516  const GObservation& obs)
517 {
518  // Get pointer on CTA observation
519  const GCTAObservation* cta = dynamic_cast<const GCTAObservation*>(&obs);
520 
521  // If pointer is not valid then throw an exception
522  if (cta == NULL) {
523  std::string cls = std::string(typeid(&obs).name());
524  std::string msg = "Invalid observation type \""+cls+"\" provided with "
525  "name \""+obs.name()+"\" (ID="+obs.id()+"). "
526  "Please specify a \"GCTAObservation\" instance as "
527  "argument.";
528  throw GException::invalid_argument(origin, msg);
529  }
530 
531  // Return reference
532  return (*cta);
533 }
534 
535 
536 /***********************************************************************//**
537  * @brief Retrieve CTA pointing from generic observation
538  *
539  * @param[in] origin Method asking for pointer retrieval.
540  * @param[in] obs Generic observation.
541  * @return Reference to CTA pointing.
542  *
543  * Extract CTA pointing from a CTA observation.
544  ***************************************************************************/
545 const GCTAPointing& gammalib::cta_pnt(const std::string& origin,
546  const GObservation& obs)
547 {
548  // Retrieve CTA observation
549  const GCTAObservation& cta = gammalib::cta_obs(origin, obs);
550 
551  // Return CTA pointing
552  return (cta.pointing());
553 }
554 
555 
556 /***********************************************************************//**
557  * @brief Retrieve CTA response from generic observation
558  *
559  * @param[in] origin Method asking for pointer retrieval.
560  * @param[in] rsp Generic response.
561  * @return Pointer to CTA response.
562  *
563  * @exception GException::invalid_argument
564  * Response @p rsp is not a CTA response.
565  *
566  * Extract CTA response from a CTA observation.
567  ***************************************************************************/
568 const GCTAResponse* gammalib::cta_rsp(const std::string& origin,
569  const GResponse& rsp)
570 {
571  // Retrieve CTA response
572  const GCTAResponse* ptr = dynamic_cast<const GCTAResponse*>(&rsp);
573 
574  // If pointer is not valid then throw an exception
575  if (ptr == NULL) {
576  std::string cls = std::string(typeid(&rsp).name());
577  std::string msg = "Invalid response type \""+cls+"\" provided. Please "
578  "specify a \"GCTAResponse\" as argument.";
579  throw GException::invalid_argument(origin, msg);
580  }
581 
582  // Return CTA response
583  return (ptr);
584 }
585 
586 
587 /***********************************************************************//**
588  * @brief Retrieve CTA IRF response from generic observation
589  *
590  * @param[in] origin Method asking for pointer retrieval.
591  * @param[in] obs Generic observation.
592  * @return Reference to CTA IRF response.
593  *
594  * @exception GException::invalid_argument
595  * Observation @p obs does not contain a CTA IRF response.
596  *
597  * Extract CTA IRF response from a CTA observation.
598  ***************************************************************************/
599 const GCTAResponseIrf& gammalib::cta_rsp_irf(const std::string& origin,
600  const GObservation& obs)
601 {
602  // Retrieve CTA observation
603  const GCTAObservation& cta = gammalib::cta_obs(origin, obs);
604 
605  // Get pointer on CTA IRF response
606  const GCTAResponseIrf* rsp = dynamic_cast<const GCTAResponseIrf*>(cta.response());
607 
608  // If pointer is not valid then throw an exception
609  if (rsp == NULL) {
610  std::string cls = std::string(typeid(cta.response()).name());
611  std::string msg = "Invalid response type \""+cls+"\" provided in "
612  "CTA observation \""+obs.name()+"\" (ID="+obs.id()+"). "
613  "Please specify a CTA observation containing an IRF "
614  "response as argument.";
615  throw GException::invalid_argument(origin, msg);
616  }
617 
618  // Return CTA IRF response
619  return (*rsp);
620 }
621 
622 
623 /***********************************************************************//**
624  * @brief Retrieve CTA cube response from generic observation
625  *
626  * @param[in] origin Method asking for pointer retrieval.
627  * @param[in] obs Generic observation.
628  * @return Reference to CTA cube response.
629  *
630  * @exception GException::invalid_argument
631  * Observation @p obs does not contain a CTA cube response.
632  *
633  * Extract CTA cube response from a CTA observation.
634  ***************************************************************************/
635 const GCTAResponseCube& gammalib::cta_rsp_cube(const std::string& origin,
636  const GObservation& obs)
637 {
638  // Retrieve CTA observation
639  const GCTAObservation& cta = gammalib::cta_obs(origin, obs);
640 
641  // Get pointer on CTA response cube
642  const GCTAResponseCube* rsp = dynamic_cast<const GCTAResponseCube*>(cta.response());
643 
644  // If pointer is not valid then throw an exception
645  if (rsp == NULL) {
646  std::string cls = std::string(typeid(cta.response()).name());
647  std::string msg = "Invalid response type \""+cls+"\" provided. Please "
648  "specify a CTA observation containing a response cube "
649  "as argument.";
650  throw GException::invalid_argument(origin, msg);
651  }
652 
653  // Return CTA response cube
654  return (*rsp);
655 }
656 
657 
658 /***********************************************************************//**
659  * @brief Retrieve CTA effective area response from generic observation
660  *
661  * @param[in] origin Method asking for pointer retrieval.
662  * @param[in] obs Generic observation.
663  * @return Reference to CTA effective area response.
664  *
665  * @exception GException::invalid_argument
666  * Observation @p obs does not contain a CTA effective area
667  * response.
668  *
669  * Extract CTA effective area response from a CTA observation.
670  ***************************************************************************/
671 const GCTAAeff& gammalib::cta_rsp_aeff(const std::string& origin,
672  const GObservation& obs)
673 {
674  // Retrieve CTA IRF response
675  const GCTAResponseIrf& rsp = gammalib::cta_rsp_irf(origin, obs);
676 
677  // Get pointer on CTA effective area
678  const GCTAAeff* aeff = rsp.aeff();
679  if (aeff == NULL) {
680  std::string msg = "Specified observation does not contain a valid "
681  "effective area response. Please specify a CTA "
682  "observation containing an effective area response.";
683  throw GException::invalid_argument(origin, msg);
684  }
685 
686  // Return CTA effective area response
687  return (*aeff);
688 }
689 
690 
691 /***********************************************************************//**
692  * @brief Retrieve CTA background response from generic observation
693  *
694  * @param[in] origin Method asking for pointer retrieval.
695  * @param[in] obs Generic observation.
696  * @return Reference to CTA background response.
697  *
698  * @exception GException::invalid_argument
699  * Observation @p obs does not contain a CTA background response.
700  *
701  * Extract CTA background response from a CTA observation.
702  ***************************************************************************/
703 const GCTABackground& gammalib::cta_rsp_bkg(const std::string& origin,
704  const GObservation& obs)
705 {
706  // Retrieve CTA IRF response
707  const GCTAResponseIrf& rsp = gammalib::cta_rsp_irf(origin, obs);
708 
709  // Get pointer on CTA background response
710  const GCTABackground* bkg = rsp.background();
711  if (bkg == NULL) {
712  std::string msg = "Specified observation does not contain a valid "
713  "background response. Please specify a CTA "
714  "observation containing a background response.";
715  throw GException::invalid_argument(origin, msg);
716  }
717 
718  // Return CTA background response
719  return (*bkg);
720 }
721 
722 
723 /***********************************************************************//**
724  * @brief Retrieve CTA event list from generic observation
725  *
726  * @param[in] origin Method asking for pointer retrieval.
727  * @param[in] obs Generic observation.
728  * @return Reference to CTA event list.
729  *
730  * @exception GException::invalid_argument
731  * Observation @p obs does not contain a CTA event list.
732  *
733  * Extract CTA event list from a CTA observation.
734  ***************************************************************************/
735 const GCTAEventList& gammalib::cta_event_list(const std::string& origin,
736  const GObservation& obs)
737 {
738  // Retrieve CTA observation
739  const GCTAObservation& cta = gammalib::cta_obs(origin, obs);
740 
741  // Get pointer on CTA events list
742  const GCTAEventList* events = dynamic_cast<const GCTAEventList*>(cta.events());
743 
744  // If pointer is not valid then throw an exception
745  if (events == NULL) {
746  std::string msg = "Specified observation does not contain a CTA event "
747  "list. Please specify a CTA observation containing "
748  "a CTA event list as argument.";
749  throw GException::invalid_argument(origin, msg);
750  }
751 
752  // Return CTA event list
753  return (*events);
754 }
755 
756 
757 /***********************************************************************//**
758  * @brief Retrieve CTA instrument direction from generic event
759  *
760  * @param[in] origin Method asking for pointer retrieval.
761  * @param[in] event Generic event.
762  * @return Reference to CTA RoI.
763  *
764  * @exception GException::invalid_argument
765  * @p event does not contain a CTA instrument direction.
766  *
767  * Extract CTA Instrument Direction from an event.
768  ***************************************************************************/
769 const GCTAInstDir& gammalib::cta_dir(const std::string& origin,
770  const GEvent& event)
771 {
772  // Get pointer on CTA instrument direction
773  const GCTAInstDir* dir = dynamic_cast<const GCTAInstDir*>(&(event.dir()));
774 
775  // If pointer is not valid then throw an exception
776  if (dir == NULL) {
777  std::string cls = std::string(typeid(&event).name());
778  std::string msg = "Invalid event type \""+cls+"\" provided. Please "
779  "specify a \"GCTAEventAtom\" or \"GCTAEventBin\" "
780  "instance as argument.";
781  throw GException::invalid_argument(origin, msg);
782  }
783 
784  // Return reference
785  return (*dir);
786 }
787 
788 
789 /***********************************************************************//**
790  * @brief Retrieve spectral component from CTA background model
791  *
792  * @param[in] model Data model.
793  * @return Pointer to spectral component.
794  *
795  * Retrieves the spectral component from a CTA background model. If the
796  * background model has no spectral component, a NULL pointer will be
797  * returned.
798  ***************************************************************************/
800 {
801  // Initialise pointer on spectral component
802  GModelSpectral* spectral = NULL;
803 
804  // Cast pointer to possible CTA background models
805  const GCTAModelBackground* cta = dynamic_cast<const GCTAModelBackground*>(&model);
806  const GCTAModelAeffBackground* aeff = dynamic_cast<const GCTAModelAeffBackground*>(&model);
807  const GCTAModelIrfBackground* irf = dynamic_cast<const GCTAModelIrfBackground*>(&model);
808  const GCTAModelCubeBackground* cube = dynamic_cast<const GCTAModelCubeBackground*>(&model);
809  const GCTAModelRadialAcceptance* rad = dynamic_cast<const GCTAModelRadialAcceptance*>(&model);
810 
811  // Extract spectral component
812  if (cta != NULL) {
813  spectral = cta->spectral();
814  }
815  else if (aeff != NULL) {
816  spectral = aeff->spectral();
817  }
818  else if (irf != NULL) {
819  spectral = irf->spectral();
820  }
821  else if (cube != NULL) {
822  spectral = cube->spectral();
823  }
824  else if (rad != NULL) {
825  spectral = rad->spectral();
826  }
827 
828  // Return pointer to spectral component
829  return spectral;
830 }
831 
832 
833 /***********************************************************************//**
834  * @brief Retrieve temporal component from CTA background model
835  *
836  * @param[in] model Data model.
837  * @return Pointer to temporal component.
838  *
839  * Retrieves the temporal component from a CTA background model. If the
840  * background model has no temporal component, a NULL pointer will be
841  * returned.
842  ***************************************************************************/
844 {
845  // Initialise pointer on temporal component
846  GModelTemporal* temporal = NULL;
847 
848  // Cast pointer to possible CTA background models
849  const GCTAModelBackground* cta = dynamic_cast<const GCTAModelBackground*>(&model);
850  const GCTAModelAeffBackground* aeff = dynamic_cast<const GCTAModelAeffBackground*>(&model);
851  const GCTAModelIrfBackground* irf = dynamic_cast<const GCTAModelIrfBackground*>(&model);
852  const GCTAModelCubeBackground* cube = dynamic_cast<const GCTAModelCubeBackground*>(&model);
853  const GCTAModelRadialAcceptance* rad = dynamic_cast<const GCTAModelRadialAcceptance*>(&model);
854 
855  // Extract temporal component
856  if (cta != NULL) {
857  temporal = cta->temporal();
858  }
859  else if (aeff != NULL) {
860  temporal = aeff->temporal();
861  }
862  else if (irf != NULL) {
863  temporal = irf->temporal();
864  }
865  else if (cube != NULL) {
866  temporal = cube->temporal();
867  }
868  else if (rad != NULL) {
869  temporal = rad->temporal();
870  }
871 
872  // Return pointer to temporal component
873  return temporal;
874 }
CTA background model base class definition.
GCTARoi read_ds_roi(const GFitsHDU &hdu)
Extract ROI from data sub-space keywords.
Definition: GCTASupport.cpp:91
std::string read_ds_gti_extname(const GFitsHDU &hdu)
Return Good Time Intervals extension name from data sub-space keywords.
CTA event list class interface definition.
void append(const double &pmin, const double &pmax)
Append phase interval.
Definition: GPhases.cpp:198
CTA Aeff background model class definition.
bool has_card(const int &cardno) const
Check existence of header card.
Definition: GFitsHDU.hpp:233
const double inv_ln2
Definition: GMath.hpp:47
const GCTAInstDir & cta_dir(const std::string &origin, const GEvent &event)
Retrieve CTA instrument direction from generic event.
CTA cube-style response function class definition.
Abstract spectral model base class.
CTA IRF background model class.
CTA IRF background model class definition.
GModelSpectral * spectral(void) const
Return spectral model component.
Abstract temporal model base class.
Abstract FITS extension base class.
Definition: GFitsHDU.hpp:51
void append(const GEnergy &emin, const GEnergy &emax)
Append energy interval.
Definition: GEbounds.cpp:302
const GCTABackground * background(void) const
Return pointer to background model.
CTA cube-style response function class.
GModelTemporal * temporal(void) const
Return temporal model component.
CTA instrument response function class definition.
const GCTAPointing & cta_pnt(const std::string &origin, const GObservation &obs)
Retrieve CTA pointing from generic observation.
Abstract interface for the event classes.
Definition: GEvent.hpp:71
CTA event list class.
const double & radius(void) const
Returns radius of region of interest in degrees.
Definition: GCTARoi.hpp:124
CTA pointing class interface definition.
virtual void response(const GResponse &rsp)
Set response function.
Radial acceptance model class interface definition.
GModelSpectral * spectral(void) const
Return spectral model component.
std::string strip_chars(const std::string &arg, const std::string &chars)
Strip leading and trailing character from string.
Definition: GTools.cpp:87
std::vector< std::string > split(const std::string &s, const std::string &sep)
Split string.
Definition: GTools.cpp:862
Interface for the CTA region of interest class.
Definition: GCTARoi.hpp:48
const GCTAResponseCube & cta_rsp_cube(const std::string &origin, const GObservation &obs)
Retrieve CTA cube response from generic observation.
Gammalib tools definition.
FITS file class.
Definition: GFits.hpp:63
#define G_READ_DS_ROI
Definition: GCTASupport.cpp:62
GPhases read_ds_phase(const GFitsHDU &hdu)
Read phase boundary data sub-space keywords.
void pointing(const GCTAPointing &pointing)
Set CTA pointing.
std::string strip_whitespace(const std::string &arg)
Strip leading and trailing whitespace from string.
Definition: GTools.cpp:73
const GModelTemporal * cta_model_temporal(const GModelData &model)
Retrieve temporal component from CTA background model.
const GCTABackground & cta_rsp_bkg(const std::string &origin, const GObservation &obs)
Retrieve CTA background response from generic observation.
void id(const std::string &id)
Set observation identifier.
Abstract FITS extension base class definition.
Definition of support function used by CTA classes.
double real(const std::string &keyname) const
Return card value as double precision.
Definition: GFitsHDU.hpp:423
CTA effective area base class definition.
Abstract temporal model base class interface definition.
GModelSpectral * spectral(void) const
Return spectral model component.
Phase intervals class interface definition.
Energy boundaries container class.
Definition: GEbounds.hpp:60
GModelTemporal * temporal(void) const
Return temporal model component.
const GCTAResponse * cta_rsp(const std::string &origin, const GResponse &rsp)
Retrieve CTA response from generic observation.
const GCTAAeff * aeff(void) const
Return pointer to effective area response.
const GCTAInstDir & centre(void) const
Returns region of interest centre.
Definition: GCTARoi.hpp:110
Abstract data model base class interface definition.
GCTAInstDir instdir(const GSkyDir &skydir) const
Get instrument direction from sky direction.
GVector log(const GVector &vector)
Computes natural logarithm of vector elements.
Definition: GVector.cpp:1184
CTA cube background model class.
CTA pointing class.
Abstract data model class.
Definition: GModelData.hpp:53
#define G_READ_DS_GTI
Definition: GCTASupport.cpp:65
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.
const GModelSpectral * cta_model_spectral(const GModelData &model)
Retrieve spectral component from CTA background model.
double resolution(const GModelSpatial *model)
Determine resolution of spatial model.
int integer(const std::string &keyname) const
Return card value as integer.
Definition: GFitsHDU.hpp:436
GModelSpectral * spectral(void) const
Return spectral model component.
const GCTAObservation & cta_obs(const std::string &origin, const GObservation &obs)
Retrieve CTA observation from generic observation.
Abstract observation base class.
void radec_deg(const double &ra, const double &dec)
Set equatorial sky direction (degrees)
Definition: GSkyDir.cpp:202
CTA response abstract base class definition.
const std::string & extname(void) const
Return extension name.
Definition: GFitsHDU.hpp:162
CTA region of interest class interface definition.
CTA observation class interface definition.
std::string gadf_hduclas4(const GFits &fits, const std::string &hduclas4)
Return extension name for GADF response table of given HDU class 4.
CTA instrument response function class.
CTA instrument response function class.
Abstract base class for the CTA effective area.
Definition: GCTAAeff.hpp:54
Abstract spectral model base class interface definition.
Background model class.
Abstract response base class definition.
void name(const std::string &name)
Set observation name.
GModelTemporal * temporal(void) const
Return temporal model component.
int iter_phi(const double &rho, const double &resolution, const int &iter_min, const int &iter_max)
Determine number of azimuthal Romberg iterations.
GModelTemporal * temporal(void) const
Return temporal model component.
CTA cube background model class interface definition.
GModelSpectral * spectral(void) const
Return spectral model component.
int size(void) const
Return number of HDUs in FITS file.
Definition: GFits.hpp:237
virtual GEvents * events(void)
Return events.
std::string string(const std::string &keyname) const
Return card value as string.
Definition: GFitsHDU.hpp:410
Energy boundaries class interface definition.
Abstract base class for the CTA background model.
Radial acceptance model class.
Phase Intervals class.
Definition: GPhases.hpp:42
CTA instrument direction class.
Definition: GCTAInstDir.hpp:59
const double twopi
Definition: GMath.hpp:36
CTA observation class.
Abstract instrument response base class.
Definition: GResponse.hpp:67
#define G_READ_DS_EBOUNDS
Definition: GCTASupport.cpp:63
GEbounds read_ds_ebounds(const GFitsHDU &hdu)
Read energy boundary data sub-space keywords.
const GCTAAeff & cta_rsp_aeff(const std::string &origin, const GObservation &obs)
Retrieve CTA effective area response from generic observation.
int iter_rho(const double &rho_max, const double &resolution, const int &iter_min, const int &iter_max)
Determine number of radial Romberg iterations.
Sky direction class.
Definition: GSkyDir.hpp:62
#define G_READ_DS_PHASE
Definition: GCTASupport.cpp:64
const GCTAResponseIrf & cta_rsp_irf(const std::string &origin, const GObservation &obs)
Retrieve CTA IRF response from generic observation.
GModelTemporal * temporal(void) const
Return temporal model component.
Mathematical function definitions.
double todouble(const std::string &arg)
Convert string into double precision value.
Definition: GTools.cpp:805
Class that handles energies in a unit independent way.
Definition: GEnergy.hpp:48
Background model class interface definition.
std::string str(const unsigned short int &value)
Convert unsigned short integer value into string.
Definition: GTools.cpp:413