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