GammaLib  2.1.0.dev
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
GCTAObservation.cpp
Go to the documentation of this file.
1 /***************************************************************************
2  * GCTAObservation.cpp - CTA Observation class *
3  * ----------------------------------------------------------------------- *
4  * copyright (C) 2010-2022 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 GCTAObservation.cpp
23  * @brief CTA observation class implementation
24  * @author Juergen Knoedlseder
25  */
26 
27 /* __ Includes ___________________________________________________________ */
28 #ifdef HAVE_CONFIG_H
29 #include <config.h>
30 #endif
31 #include <typeinfo>
32 #include "GObservationRegistry.hpp"
33 #include "GException.hpp"
34 #include "GFits.hpp"
35 #include "GFilename.hpp"
36 #include "GGti.hpp"
37 #include "GTools.hpp"
38 #include "GIntegral.hpp"
39 #include "GCTASupport.hpp"
40 #include "GCTAObservation.hpp"
41 #include "GCTAResponseIrf.hpp"
42 #include "GCTAResponseCube.hpp"
43 #include "GCTAEventList.hpp"
44 #include "GCTAEventCube.hpp"
45 #include "GCTARoi.hpp"
46 
47 /* __ Globals ____________________________________________________________ */
48 const GCTAObservation g_obs_cta_seed(true, "CTA");
49 const GCTAObservation g_obs_hess_seed(true, "HESS");
50 const GCTAObservation g_obs_magic_seed(true, "MAGIC");
51 const GCTAObservation g_obs_veritas_seed(true, "VERITAS");
52 const GCTAObservation g_obs_astri_seed(true, "ASTRI");
53 const GCTAObservation g_obs_fact_seed(true, "FACT");
54 const GObservationRegistry g_obs_cta_registry(&g_obs_cta_seed);
55 const GObservationRegistry g_obs_hess_registry(&g_obs_hess_seed);
56 const GObservationRegistry g_obs_magic_registry(&g_obs_magic_seed);
57 const GObservationRegistry g_obs_veritas_registry(&g_obs_veritas_seed);
58 const GObservationRegistry g_obs_astri_registry(&g_obs_astri_seed);
59 const GObservationRegistry g_obs_fact_registry(&g_obs_fact_seed);
60 
61 /* __ Method name definitions ____________________________________________ */
62 #define G_RESPONSE_SET "GCTAObservation::response(GResponse&)"
63 #define G_RESPONSE_GET "GCTAObservation::response()"
64 #define G_ROI "GCTAObservation::roi()"
65 #define G_GTI "GCTAObservation::gti()"
66 #define G_EBOUNDS "GCTAObservation::ebounds()"
67 #define G_READ "GCTAObservation::read(GXmlElement&)"
68 #define G_WRITE "GCTAObservation::write(GXmlElement&)"
69 #define G_LOAD1 "GCTAObservation::load(GFilename&, GFilename&, "\
70  "GFilename&, GFilename&)"
71 #define G_LOAD2 "GCTAObservation::load(GFilename&, GFilename&, "\
72  "GFilename&, GFilename&, GFilename&)"
73 
74 /* __ Macros _____________________________________________________________ */
75 
76 /* __ Coding definitions _________________________________________________ */
77 
78 /* __ Debug definitions __________________________________________________ */
79 
80 
81 /*==========================================================================
82  = =
83  = Constructors/destructors =
84  = =
85  ==========================================================================*/
86 
87 /***********************************************************************//**
88  * @brief Void constructor
89  *
90  * Constructs an empty CTA observation.
91  ***************************************************************************/
93 {
94  // Initialise members
95  init_members();
96 
97  // Return
98  return;
99 }
100 
101 
102 /***********************************************************************//**
103  * @brief Instrument constructor
104  *
105  * @param[in] dummy Dummy flag.
106  * @param[in] instrument Instrument name.
107  *
108  * Constructs an empty CTA observation for a given instrument.
109  *
110  * This method is specifically used allocation of global class instances for
111  * observation registry. By specifying explicit instrument names it is
112  * possible to use the "CTA" module are for other Imaging Air Cherenkov
113  * Telescopes. So far, the following instrument codes are supported:
114  * "CTA", "HESS", "VERITAS", "MAGIC", "ASTRI".
115  ***************************************************************************/
117  const std::string& instrument) :
118  GObservation()
119 {
120  // Initialise members
121  init_members();
122 
123  // Set instrument
125 
126  // Return
127  return;
128 }
129 
130 
131 /***********************************************************************//**
132  * @brief Unbinned or binned analysis constructor
133  *
134  * @param[in] filename Event list or counts cube FITS file name.
135  *
136  * Constructs a CTA observation from either an event list of a counts cube.
137  ***************************************************************************/
139 {
140  // Initialise members
141  init_members();
142 
143  // Load data
144  load(filename);
145 
146  // Return
147  return;
148 }
149 
150 
151 /***********************************************************************//**
152  * @brief Stacked analysis constructor
153  *
154  * @param[in] cntcube Counts cube file name.
155  * @param[in] expcube Exposure cube file name.
156  * @param[in] psfcube Point spread function cube file name.
157  * @param[in] bkgcube Background cube file name.
158  *
159  * Constructs a CTA observation from a counts cube, an exposure cube, a point
160  * spread function cube, and a background cube.
161  ***************************************************************************/
163  const GFilename& expcube,
164  const GFilename& psfcube,
165  const GFilename& bkgcube) : GObservation()
166 {
167  // Initialise members
168  init_members();
169 
170  // Load data
171  load(cntcube, expcube, psfcube, bkgcube);
172 
173  // Return
174  return;
175 }
176 
177 
178 /***********************************************************************//**
179  * @brief Stacked analysis with energy dispersion constructor
180  *
181  * @param[in] cntcube Counts cube file name.
182  * @param[in] expcube Exposure cube file name.
183  * @param[in] psfcube Point spread function cube file name.
184  * @param[in] edispcube Energy dispersion cube file name.
185  * @param[in] bkgcube Background cube file name.
186  *
187  * Constructs a CTA observation from a counts cube, an exposure cube, a point
188  * spread function cube, an energy dispersion cube, and a background cube.
189  ***************************************************************************/
191  const GFilename& expcube,
192  const GFilename& psfcube,
193  const GFilename& edispcube,
194  const GFilename& bkgcube) : GObservation()
195 {
196  // Initialise members
197  init_members();
198 
199  // Load data
200  load(cntcube, expcube, psfcube, edispcube, bkgcube);
201 
202  // Return
203  return;
204 }
205 
206 
207 /***********************************************************************//**
208  * @brief Copy constructor
209  *
210  * @param[in] obs CTA observation.
211  *
212  * Constructs a CTA observation by copying an existing CTA observation.
213  ***************************************************************************/
215 {
216  // Initialise members
217  init_members();
218 
219  // Copy members
220  copy_members(obs);
221 
222  // Return
223  return;
224 }
225 
226 
227 /***********************************************************************//**
228  * @brief Destructor
229  *
230  * Destructs CTA observation.
231  ***************************************************************************/
233 {
234  // Free members
235  free_members();
236 
237  // Return
238  return;
239 }
240 
241 
242 /*==========================================================================
243  = =
244  = Operators =
245  = =
246  ==========================================================================*/
247 
248 /***********************************************************************//**
249  * @brief Assignment operator
250  *
251  * @param[in] obs CTA observation.
252  *
253  * Assign CTA observation to this object.
254  ***************************************************************************/
256 {
257  // Execute only if object is not identical
258  if (this != &obs) {
259 
260  // Copy base class members
261  this->GObservation::operator=(obs);
262 
263  // Free members
264  free_members();
265 
266  // Initialise members
267  init_members();
268 
269  // Copy members
270  copy_members(obs);
271 
272  } // endif: object was not identical
273 
274  // Return this object
275  return *this;
276 }
277 
278 
279 /*==========================================================================
280  = =
281  = Public methods =
282  = =
283  ==========================================================================*/
284 
285 /***********************************************************************//**
286  * @brief Clear CTA observation
287  *
288  * Clear CTA observation.
289  ***************************************************************************/
291 {
292  // Free members
293  free_members();
295 
296  // Initialise members
298  init_members();
299 
300  // Return
301  return;
302 }
303 
304 
305 /***********************************************************************//**
306  * @brief Clone CTA observation
307  *
308  * @return Pointer to deep copy of CTA observation.
309  *
310  * Returns a pointer to a deep copy of a CTA observation.
311  ***************************************************************************/
313 {
314  return new GCTAObservation(*this);
315 }
316 
317 
318 /***********************************************************************//**
319  * @brief Set response function
320  *
321  * @param[in] rsp Response function.
322  *
323  * Sets the response function for the observation.
324  ***************************************************************************/
326 {
327  // Retrieve CTA response pointer
328  const GCTAResponse* ptr = gammalib::cta_rsp(G_RESPONSE_SET, rsp);
329 
330  // Free existing response only if it differs from current response. This
331  // prevents unintential deallocation of the argument
332  if ((m_response != NULL) && (m_response != &rsp)) {
333  delete m_response;
334  }
335 
336  // Clone response function
337  m_response = ptr->clone();
338 
339  // Return
340  return;
341 }
342 
343 
344 /***********************************************************************//**
345  * @brief Return pointer to CTA response function
346  *
347  * @return Pointer to CTA response function.
348  *
349  * @exception GException::invalid_value
350  * No valid response found in CTA observation.
351  *
352  * Returns a pointer to the CTA response function. An exception is thrown if
353  * the pointer is not valid, hence the user does not need to verify the
354  * validity of the pointer.
355  ***************************************************************************/
357 {
358  // Throw an exception if the response pointer is not valid
359  if (m_response == NULL) {
360  std::string msg = "No valid response function found in CTA "
361  "observation.\n";
363  }
364 
365  // Return pointer
366  return m_response;
367 }
368 
369 
370 /***********************************************************************//**
371  * @brief Set CTA response function
372  *
373  * @param[in] rspname Response name.
374  * @param[in] caldb Calibration database.
375  *
376  * Sets the CTA response function by specifying a response name and a
377  * calibration database. This method also loads the response function so that
378  * it is available for data analysis.
379  ***************************************************************************/
380 void GCTAObservation::response(const std::string& rspname, const GCaldb& caldb)
381 {
382  // Free response
383  if (m_response != NULL) delete m_response;
384  m_response = NULL;
385 
386  // Allocate fresh response function
387  GCTAResponseIrf* rsp = new GCTAResponseIrf;
388 
389  // Set calibration database
390  rsp->caldb(caldb);
391 
392  // Load instrument response function
393  rsp->load(rspname);
394 
395  // Store pointer
396  m_response = rsp;
397 
398  // Return
399  return;
400 }
401 
402 
403 /***********************************************************************//**
404  * @brief Set CTA response function
405  *
406  * @param[in] expcube Exposure cube.
407  * @param[in] psfcube Psf cube.
408  * @param[in] bkgcube Background cube.
409  *
410  * Sets the CTA response function fur cube analysis by specifying the
411  * exposure cube, the Psf cube and the background cube. The method also
412  * copies over the ontime, the livetime and the deadtime correction factor
413  * from the exposure cube.
414  ***************************************************************************/
416  const GCTACubePsf& psfcube,
417  const GCTACubeBackground& bkgcube)
418 {
419  // Free response
420  if (m_response != NULL) delete m_response;
421  m_response = NULL;
422 
423  // Allocate fresh response function
424  GCTAResponseCube* rsp = new GCTAResponseCube(expcube, psfcube, bkgcube);
425 
426  // Store pointer
427  m_response = rsp;
428 
429  // Copy over time information from exposure cube
430  ontime(expcube.ontime());
431  livetime(expcube.livetime());
432  deadc(expcube.deadc());
433 
434  // Return
435  return;
436 }
437 
438 
439 /***********************************************************************//**
440  * @brief Set CTA response function
441  *
442  * @param[in] expcube Exposure cube.
443  * @param[in] psfcube Psf cube.
444  * @param[in] edispcube Edisp cube.
445  * @param[in] bkgcube Background cube.
446  *
447  * Sets the CTA response function fur cube analysis by specifying the
448  * exposure cube, the Psf cube, the exposure cube and the background cube.
449  * The method also copies over the ontime, the livetime and the deadtime
450  * correction factor from the exposure cube.
451  ***************************************************************************/
453  const GCTACubePsf& psfcube,
454  const GCTACubeEdisp& edispcube,
455  const GCTACubeBackground& bkgcube)
456 {
457  // Free response
458  if (m_response != NULL) delete m_response;
459  m_response = NULL;
460 
461  // Allocate fresh response function
462  GCTAResponseCube* rsp = new GCTAResponseCube(expcube,
463  psfcube,
464  edispcube,
465  bkgcube);
466 
467  // Store pointer
468  m_response = rsp;
469 
470  // Copy over time information from exposure cube
471  ontime(expcube.ontime());
472  livetime(expcube.livetime());
473  deadc(expcube.deadc());
474 
475  // Return
476  return;
477 }
478 
479 
480 /***********************************************************************//**
481  * @brief Get Region of Interest
482  *
483  * @return Region of Interest.
484  *
485  * @exception GException::invalid_value
486  * Observation does not contain events.
487  * Observation does not contain an event list.
488  *
489  * Extracts the Region of Interest from the event list. An exception is
490  * thrown if the observation does not contain an event list.
491  ***************************************************************************/
493 {
494  // Throw an exception if no events exist
495  if (m_events == NULL) {
496  std::string msg = "Region of Interest is not defined since the "
497  "observation does not contain events.";
498  throw GException::invalid_value(G_ROI, msg);
499  }
500 
501  // Get pointer to event list. Throw an exception if no event list is found
502  const GCTAEventList* list = dynamic_cast<const GCTAEventList*>(m_events);
503  if (list == NULL) {
504  std::string msg = "Region of Interest is not defined since the "
505  "observation does not contain an event list.";
506  throw GException::invalid_value(G_ROI, msg);
507  }
508 
509  // Return ROI
510  return (list->roi());
511 }
512 
513 
514 /***********************************************************************//**
515  * @brief Get Good Time Intervals
516  *
517  * @return Good Time Intervals.
518  *
519  * @exception GException::invalid_value
520  * Observation does not contain events.
521  *
522  * Extracts the Good Time Intervals from the events. An exception is thrown
523  * if the observation does not contain events.
524  ***************************************************************************/
526 {
527  // Throw an exception if no events exist
528  if (m_events == NULL) {
529  std::string msg = "Good Time Intervals are not defined since the "
530  "observation does not contain events.";
531  throw GException::invalid_value(G_GTI, msg);
532  }
533 
534  // Return GTI
535  return (m_events->gti());
536 }
537 
538 
539 /***********************************************************************//**
540  * @brief Get energy boundaries
541  *
542  * @return Energy boundaries.
543  *
544  * @exception GException::invalid_value
545  * Observation does not contain events.
546  *
547  * Extract the energy boundaries from the events. An exception is thrown if
548  * the observation does not contain events.
549  ***************************************************************************/
551 {
552  // Throw an exception if no events exist
553  if (m_events == NULL) {
554  std::string msg = "Energy boundaries are not defined since the "
555  "observation does not contain events.";
557  }
558 
559  // Return energy boundaries
560  return (m_events->ebounds());
561 }
562 
563 
564 /***********************************************************************//**
565  * @brief Read observation from XML element
566  *
567  * @param[in] xml XML element.
568  *
569  * @exception GException::invalid_value
570  * Invalid number of parameters found in XML element.
571  * @exception GException::xml_invalid_parnames
572  * Invalid parameter names found in XML element.
573  * @exception GException::invalid_value
574  * Invalid statistic attribute encountered
575  *
576  * Reads information for a CTA observation from an XML element. This method
577  * handles two variants: a first where an event list of counts cube is
578  * given and a second where the observation definition information is
579  * provided by parameter tags.
580  *
581  * The XML format for an event list is
582  *
583  * <observation name="..." id="..." instrument="...">
584  * <parameter name="EventList" file="..."/>
585  * ...
586  * </observation>
587  *
588  * and for a counts cube is
589  *
590  * <observation name="..." id="..." instrument="...">
591  * <parameter name="CountsCube" file="..."/>
592  * ...
593  * </observation>
594  *
595  * The second variant without event information has the following format
596  *
597  * <observation name="..." id="..." instrument="...">
598  * <parameter name="Pointing" ra="..." dec="..."/>
599  * <parameter name="Deadtime" deadc="..."/>
600  * ...
601  * </observation>
602  *
603  * In addition, calibration information can be specified using the format
604  *
605  * <observation name="..." id="..." instrument="...">
606  * ...
607  * <parameter name="Calibration" database="..." response="..."/>
608  * </observation>
609  *
610  * If even more control is required over the response, individual file names
611  * can be specified using
612  *
613  * <observation name="..." id="..." instrument="...">
614  * ...
615  * <parameter name="EffectiveArea" file="..."/>
616  * <parameter name="PointSpreadFunction" file="..."/>
617  * <parameter name="EnergyDispersion" file="..."/>
618  * <parameter name="Background" file="..."/>
619  * </observation>
620  *
621  * In case that a @a CountsCube is handled, the stacked response can also be
622  * provided in the format
623  *
624  * <observation name="..." id="..." instrument="...">
625  * ...
626  * <parameter name="ExposureCube" file="..."/>
627  * <parameter name="PsfCube" file="..."/>
628  * <parameter name="EdispCube" file="..."/>
629  * <parameter name="BkgCube" file="..."/>
630  * </observation>
631  *
632  * Note that the @p EdispCube is an optional parameter.
633  *
634  * Optional user energy thresholds can be specified by adding the @a emin
635  * and @a emax attributes to the @p observation tag. Also the statistic
636  * used for maximum likelihood fitting can be specified:
637  *
638  * <observation name="..." id="..." instrument="..." emin="..." emax="..." statistic="...">
639  *
640  * The method does not load the events into memory but stores the file name
641  * of the event file. The events are only loaded when required. This reduces
642  * the memory needs for an CTA observation object and allows for loading
643  * of event information upon need.
644  ***************************************************************************/
646 {
647  // Clear observation
648  clear();
649 
650  // Extract instrument name
651  m_instrument = xml.attribute("instrument");
652 
653  // Read in user defined energy boundaries of this observation
654  if (xml.attribute("emin") != "") {
656  }
657  if (xml.attribute("emax") != "") {
659  }
660 
661  // Read in user defined statistic for this observation
662  if (xml.attribute("statistic") != "") {
663 
664  // Extract statistic value
665  std::string statistic = gammalib::toupper(xml.attribute("statistic"));
666 
667  // If statistic is not POISSON, CSTAT, GAUSSIAN or CHI2 than throw
668  // an exception
669  if ((statistic != "POISSON") &&
670  (statistic != "CSTAT") &&
671  (statistic != "GAUSSIAN") &&
672  (statistic != "CHI2")) {
673  std::string msg = "Invalid statistic \""+statistic+"\" encountered "
674  "in observation definition XML file for "
675  "\""+m_instrument+"\" observation with identifier "
676  "\""+xml.attribute("id")+"\". Only \"POISSON\", "
677  "\"CSTAT\", \"GAUSSIAN\" or \"CHI2\" are "
678  "supported.";
679  throw GException::invalid_value(G_READ, msg);
680  }
681 
682  // Save statistic value
683  this->statistic(xml.attribute("statistic"));
684 
685  }
686 
687  // Determine number of parameter nodes in XML element
688  int npars = xml.elements("parameter");
689 
690  // Verify that XML element has at least 1 parameter
691  if (xml.elements() < 1 || npars < 1) {
692  std::string msg = "CTA observation contains "+gammalib::str(npars)+
693  " parameters but at least one parameter is required. "
694  "Please verify the XML format.";
695  throw GException::invalid_value(G_READ, msg);
696  }
697 
698  // First try to extract observation parameters for an event file
699  int n_evfile = 0;
700  for (int i = 0; i < npars; ++i) {
701 
702  // Get parameter element
703  const GXmlElement* par = xml.element("parameter", i);
704 
705  // Handle EventList or CountsCube
706  if ((par->attribute("name") == "EventList") ||
707  (par->attribute("name") == "CountsCube")) {
708 
709  // Read eventlist file name
710  std::string filename = gammalib::xml_file_expand(xml,
711  par->attribute("file"));
712 
713  // Load events
714  load(filename);
715 
716  // Increment parameter counter
717  n_evfile++;
718  }
719 
720  // Read Background filename (needed by GCTAModelCubeBackground)
721  else if (par->attribute("name") == "Background") {
722 
723  // Read background file name
724  m_bgdfile = gammalib::xml_file_expand(xml, par->attribute("file"));
725 
726  }
727 
728  // Read Off regions
729  else if (par->attribute("name") == "OffRegions") {
730 
731  // Read off regions file name
732  std::string filename = gammalib::xml_file_expand(xml,
733  par->attribute("file"));
734 
735  // Load off regions
736  m_off_regions = GSkyRegions(filename);
737 
738  }
739 
740  } // endfor: looped over observation parameters
741 
742  // Analyse parameters
743  bool has_evfile = (n_evfile > 0);
744 
745  // If we have an event file then verify that all required parameters
746  // were found
747  if (has_evfile) {
748  gammalib::xml_check_par(G_READ, "EventList\" or \"CountsCube", n_evfile);
749  }
750 
751  // ... otherwise extract information from observation definition parameters
752  else {
753 
754  // Initialise event list definition
756  GGti gti;
757  GCTARoi roi;
758  bool has_ebounds = false;
759  bool has_gti = false;
760  bool has_roi = false;
761 
762  // Extract parameters (do nothing if parameter does not exist)
763  if (gammalib::xml_has_par(xml, "Pointing")) {
764  m_pointing.read(xml);
765  }
766  if (gammalib::xml_has_par(xml, "EnergyBoundaries")) {
767  ebounds.read(xml);
768  has_ebounds = true;
769  }
770  if (gammalib::xml_has_par(xml, "GoodTimeIntervals")) {
771  gti.read(xml);
772  has_gti = true;
773  }
774  if (gammalib::xml_has_par(xml, "RegionOfInterest")) {
775  roi.read(xml);
776  has_roi = true;
777  }
778  if (gammalib::xml_has_par(xml, "Deadtime")) {
779  const GXmlElement* par = gammalib::xml_get_par(G_READ, xml, "Deadtime");
780  m_deadc = gammalib::todouble(par->attribute("deadc"));
781  }
782  else {
783  m_deadc = 1.0;
784  }
785 
786  // If we have at least one event list definition then create event
787  // list and attach it to the observation
788  if (has_ebounds || has_gti || has_roi) {
789 
790  // Create an empty event list
792 
793  // Handle energy boundaries
794  if (has_ebounds) {
795  events.ebounds(ebounds);
796  }
797 
798  // Handle GTIs
799  if (has_gti) {
800  events.gti(gti);
801  ontime(gti.ontime());
802  livetime(gti.ontime() * m_deadc);
803  }
804 
805  // Handle ROIs
806  if (has_roi) {
807  events.roi(roi);
808  }
809 
810  // Attach event list
811  this->events(events);
812 
813  } // endif: handled event list information
814 
815  } // endelse: extracted information from observation definition
816 
817  // Determine response type as function of the information that is
818  // provided in the XML file. The response type is determined from the
819  // parameter names. An exception is thrown if the response type cannot
820  // be unambigously determined
821  int response_type = 0;
822  for (int i = 0; i < npars; ++i) {
823 
824  // Get parameter element
825  const GXmlElement* par = xml.element("parameter", i);
826 
827  // Check for response type 1 (GCTAResponseIrf)
828  if ((par->attribute("name") == "EffectiveArea") ||
829  (par->attribute("name") == "ARF") ||
830  (par->attribute("name") == "PointSpreadFunction") ||
831  (par->attribute("name") == "PSF") ||
832  (par->attribute("name") == "EnergyDispersion") ||
833  (par->attribute("name") == "RMF") ||
834  (par->attribute("name") == "Background") ||
835  (par->attribute("name") == "Calibration")) {
836  if (response_type == 2) {
837  std::string msg = "Response cube parameters expected but the "
838  "IRF response parameter \""+
839  par->attribute("name")+"\" was encountered. "
840  "Please verify the XML format.";
841  throw GException::invalid_value(G_READ, msg);
842  }
843  response_type = 1;
844  }
845 
846  // Check for response type 2 (GCTAResponseCube)
847  else if ((par->attribute("name") == "ExposureCube") ||
848  (par->attribute("name") == "PsfCube") ||
849  (par->attribute("name") == "BkgCube")) {
850  if (response_type == 1) {
851  std::string msg = "IRF response parameters expected but the "
852  "response cube parameter \""+
853  par->attribute("name")+"\" was encountered. "
854  "Please verify the XML format.";
855  throw GException::invalid_value(G_READ, msg);
856  }
857  response_type = 2;
858  }
859 
860  } // endfor: looped over all parameters
861 
862  // Allocate response
863  switch (response_type) {
864  case 1:
866  break;
867  case 2:
869  break;
870  default:
871  break;
872  }
873 
874  // Extract response information
875  if (m_response != NULL) {
876  m_response->read(xml);
877  }
878 
879  // Return
880  return;
881 }
882 
883 
884 /***********************************************************************//**
885  * @brief Write observation into XML element
886  *
887  * @param[in] xml XML element.
888  *
889  * @exception GException::invalid_value
890  * No valid events found in observation.
891  *
892  * Writes information for a CTA observation into an XML element.
893  *
894  * Dependent on the content of the CTA observation, different XML formats
895  * will be written. If the CTA observation contains an event list that
896  * has been loaded from disk, the method will write the file name of the
897  * event list using the format
898  *
899  * <observation name="..." id="..." instrument="..." statistic="...">
900  * <parameter name="EventList" file="..."/>
901  * ...
902  * </observation>
903  *
904  * If the CTA observation contains a counts cube that has been loaded from
905  * disk, the method will write the file name of the counts cube using the
906  * format
907  *
908  * <observation name="..." id="..." instrument="..." statistic="...">
909  * <parameter name="CountsCube" file="..."/>
910  * ...
911  * </observation>
912  *
913  * In all other cases the method will only write the metadata information
914  * for the CTA observation using the format
915  *
916  * <observation name="..." id="..." instrument="..." statistic="...">
917  * <parameter name="Pointing" ra="..." dec="..."/>
918  * <parameter name="Deadtime" deadc="..."/>
919  * ...
920  * </observation>
921  *
922  * In case that response information is available and that the response
923  * information has been either loaded or saved to disk, the method will
924  * write the reponse information into the XML file. If the response for
925  * an unbinned or binned observation has been loaded from the calibration
926  * database, the format
927  *
928  * <observation name="..." id="..." instrument="..." statistic="...">
929  * ...
930  * <parameter name="Calibration" database="..." response="..."/>
931  * </observation>
932  *
933  * will be written. If the response has been loaded from individual files,
934  * the format
935  *
936  * <observation name="..." id="..." instrument="..." statistic="...">
937  * ...
938  * <parameter name="EffectiveArea" file="..."/>
939  * <parameter name="PointSpreadFunction" file="..."/>
940  * <parameter name="EnergyDispersion" file="..."/>
941  * <parameter name="Background" file="..."/>
942  * </observation>
943  *
944  * will be written. If the observation contains information about the
945  * response to a stacked analysis, the format
946  *
947  * <observation name="..." id="..." instrument="..." statistic="...">
948  * ...
949  * <parameter name="ExposureCube" file="..."/>
950  * <parameter name="PsfCube" file="..."/>
951  * <parameter name="EdispCube" file="..."/>
952  * <parameter name="BkgCube" file="..."/>
953  * </observation>
954  *
955  * will be written. Note that the @p EdispCube parameter will only be
956  * written in case that an energy dispersion cube had been defined (energy
957  * dispersion is optional).
958  *
959  * If user energy thresholds are defined (i.e. threshold values are >0) the
960  * additional @a emin and @a emax attributes will be written to the
961  * @p observation tag:
962  *
963  * <observation name="..." id="..." instrument="..." statistic="..." emin="..." emax="...">
964  *
965  ***************************************************************************/
967 {
968  // Throw an exception if eventtype() is neither "EventList" nor
969  std::string evttype = eventtype();
970  if ((evttype != "EventList") && (evttype != "CountsCube")) {
971  std::string msg;
972  if (evttype.empty()) {
973  msg = "The observation does not contain any events, hence "
974  "it cannot be written into an XML file.";
975  }
976  else {
977  msg = "The observation contains an unknown event type \""+
978  evttype+"\". The event type needs to be either "
979  "\"EventList\" or \"CountsCube\".";
980  }
982  }
983 
984  // Add user energy threshold attributes (if set)
985  if (m_lo_user_thres > 0.0) {
987  }
988  if (m_hi_user_thres > 0.0) {
990  }
991 
992  // Add user defined statistic attributes
993  if (statistic() != "") {
994  xml.attribute("statistic", statistic());
995  }
996 
997  // If there is an event filename then write the event information to the
998  // XML file ...
999  if (!m_eventfile.is_empty()) {
1000 
1001  // Write event file name
1002  GXmlElement* par = gammalib::xml_need_par(G_WRITE, xml, evttype);
1003  par->attribute("file", gammalib::xml_file_reduce(xml, m_eventfile.url()));
1004 
1005  }
1006 
1007  // ... otherwise write the observation definition information
1008  else {
1009 
1010  // Write pointing, energy bounds, GTIs and ROI
1011  m_pointing.write(xml);
1012  events()->ebounds().write(xml);
1013  events()->gti().write(xml);
1014 
1015  // Write ROI
1016  const GCTAEventList* list = dynamic_cast<const GCTAEventList*>(events());
1017  if (list != NULL) {
1018  list->roi().write(xml);
1019  }
1020 
1021  // Write deadtime correction factor
1022  GXmlElement* par = gammalib::xml_need_par(G_WRITE, xml, "Deadtime");
1023  par->attribute("deadc", gammalib::str(m_deadc));
1024 
1025  }
1026 
1027  // Write response information
1028  if (m_response != NULL) {
1029  m_response->write(xml);
1030  }
1031 
1032  // Write Off regions filename
1033  if (!m_off_regions.filename().is_empty()) {
1034  GXmlElement* par = gammalib::xml_need_par(G_WRITE, xml, "OffRegions");
1035  par->attribute("file", gammalib::xml_file_reduce(xml,
1036  m_off_regions.filename().url()));
1037  }
1038 
1039  // Return
1040  return;
1041 }
1042 
1043 
1044 /***********************************************************************//**
1045  * @brief Read data from FITS file
1046  *
1047  * @param[in] fits FITS file.
1048  *
1049  * Reads event data from a FITS file and sets the observation attributes.
1050  *
1051  * The method automatically switches between an event list and a counts cube,
1052  * depending of the information provided in the FITS file. If an extension
1053  * name is specified, the method checks whether the extension exists and
1054  * loads the extension as event list. Otherwise, it checks whether the file
1055  * contains an `EVENTS` extension and loads the extension as event list.
1056  * If none of the above are satistified, the method loads a counts cube.
1057  ***************************************************************************/
1058 void GCTAObservation::read(const GFits& fits)
1059 {
1060  // Delete any existing event container (do not call clear() as we do not
1061  // want to delete the response function)
1062  if (m_events != NULL) delete m_events;
1063  m_events = NULL;
1064 
1065  // Retrieve file name from FITS file
1066  GFilename filename(fits.filename());
1067 
1068  // Get extension name
1069  std::string extname = filename.extname(gammalib::extname_cta_events);
1070 
1071  // If FITS file contains an EVENTS extension we have an unbinned
1072  // observation ...
1073  if (fits.contains(extname)) {
1074 
1075  // Allocate event list
1077 
1078  // Assign event list as the observation's event container
1079  m_events = events;
1080 
1081  // Read event list
1082  events->read(fits);
1083 
1084  // Read observation attributes from EVENTS extension
1085  const GFitsHDU& hdu = *fits.at(extname);
1086  read_attributes(hdu);
1087 
1088  }
1089 
1090  // ... otherwise we have a binned observation
1091  else {
1092 
1093  // Allocate event cube
1095 
1096  // Assign event cube as the observation's event container
1097  m_events = events;
1098 
1099  // Read event cube
1100  events->read(fits);
1101 
1102  // Read observation attributes from primary extension
1103  const GFitsHDU& hdu = *fits.at(0);
1104  read_attributes(hdu);
1105 
1106  }
1107 
1108  // Return
1109  return;
1110 }
1111 
1112 
1113 /***********************************************************************//**
1114  * @brief Write observation into FITS file.
1115  *
1116  * @param[in] fits FITS file.
1117  * @param[in] evtname Events FITS extension name.
1118  * @param[in] gtiname Good Time Intervals FITS extension name.
1119  *
1120  * Writes the observation into a FITS file.
1121  *
1122  * If the observation contains an event list, the event list and Good Time
1123  * Intervals will be written to the FITS file. The FITS extension name for
1124  * the event list and the Good Time Intervals can be optionally specified
1125  * thru the @p evtname and @p gtiname arguments.
1126  *
1127  * If the observation contains an event cube, the event cube will be written
1128  * into the primary extension of the FITS file.
1129  *
1130  * This method does nothing if no events are in the CTA observation.
1131  ***************************************************************************/
1132 void GCTAObservation::write(GFits& fits, const std::string& evtname,
1133  const std::string& gtiname) const
1134 {
1135  // Get pointers on event list
1136  const GCTAEventList* list = dynamic_cast<const GCTAEventList*>(events());
1137  const GCTAEventCube* cube = dynamic_cast<const GCTAEventCube*>(events());
1138 
1139  // Case A: Observation contains an event list
1140  if (list != NULL) {
1141 
1142  // Remove HDUs if they exist already
1143  if (fits.contains(evtname)) {
1144  fits.remove(evtname);
1145  }
1146  if (fits.contains(gtiname)) {
1147  fits.remove(gtiname);
1148  }
1149 
1150  // Write event list and Good Time Intervals into FITS file
1151  list->write(fits, evtname, gtiname);
1152 
1153  // Get reference to events extension
1154  GFitsHDU& hdu = *fits.at(evtname);
1155 
1156  // Write observation attributes to events extension
1157  write_attributes(hdu);
1158 
1159  } // endif: observation contained an event list
1160 
1161  // Case B: Observation contains an event cube
1162  else if (cube != NULL) {
1163 
1164  // Write events cube into FITS file. This method also writes
1165  // the energy boundaries and the GTI as they are also part
1166  // of the event cube.
1167  cube->write(fits);
1168 
1169  // Write observation attributes into primary header
1170  GFitsHDU& hdu = *fits.at(0);
1171  write_attributes(hdu);
1172 
1173  } // endelse: observation contained an event cube
1174 
1175  // Return
1176  return;
1177 }
1178 
1179 
1180 /***********************************************************************//**
1181  * @brief Load unbinned or binned analysis data
1182  *
1183  * @param[in] filename Event list or counts cube FITS file name.
1184  *
1185  * Loads either an event list or a counts cube from a FITS file.
1186  ***************************************************************************/
1187 void GCTAObservation::load(const GFilename& filename)
1188 {
1189  #pragma omp critical(GCTAObservation_load)
1190  {
1191  // Store event filename
1192  m_eventfile = filename;
1193 
1194  // Open FITS file
1195  GFits fits(filename);
1196 
1197  // Read data
1198  read(fits);
1199 
1200  // Close FITS file
1201  fits.close();
1202  }
1203 
1204  // Return
1205  return;
1206 }
1207 
1208 
1209 /***********************************************************************//**
1210  * @brief Load stacked analysis data
1211  *
1212  * @param[in] cntcube Counts cube file name.
1213  * @param[in] expcube Exposure cube file name.
1214  * @param[in] psfcube Point spread function cube file name.
1215  * @param[in] bkgcube Background cube file name.
1216  *
1217  * Loads a counts cube, an exposure cube, a point spread function cube, and a
1218  * background cube for stacked analysis.
1219  ***************************************************************************/
1220 void GCTAObservation::load(const GFilename& cntcube,
1221  const GFilename& expcube,
1222  const GFilename& psfcube,
1223  const GFilename& bkgcube) {
1224 
1225  // Load counts cube FITS file
1226  load(cntcube);
1227 
1228  // Check whether we have an event cube
1229  if (dynamic_cast<const GCTAEventCube*>(events()) == NULL) {
1230  std::string msg = "Specified file \""+cntcube.url()+"\" is not a CTA "
1231  "counts cube. Please provide a counts cube.";
1233  }
1234 
1235  // Load exposure cube
1236  GCTACubeExposure exposure(expcube);
1237 
1238  // Load point spread function cube
1239  GCTACubePsf psf(psfcube);
1240 
1241  // Load background cube
1242  GCTACubeBackground background(bkgcube);
1243 
1244  // Set exposure cube, point spread function cube, and background cube
1245  // as response for this observation
1246  response(exposure, psf, background);
1247 
1248  // Return
1249  return;
1250 }
1251 
1252 
1253 /***********************************************************************//**
1254  * @brief Load stacked analysis data with energy dispersion
1255  *
1256  * @param[in] cntcube Counts cube file name.
1257  * @param[in] expcube Exposure cube file name.
1258  * @param[in] psfcube Point spread function cube file name.
1259  * @param[in] edispcube Energy dispersion cube file name.
1260  * @param[in] bkgcube Background cube file name.
1261  *
1262  * Loads a counts cube, an exposure cube, a point spread function cube, a
1263  * energy dispersion cube and a background cube for stacked analysis.
1264  ***************************************************************************/
1265 void GCTAObservation::load(const GFilename& cntcube,
1266  const GFilename& expcube,
1267  const GFilename& psfcube,
1268  const GFilename& edispcube,
1269  const GFilename& bkgcube) {
1270 
1271  // Load counts cube FITS file
1272  load(cntcube);
1273 
1274  // Check whether we have an event cube
1275  if (dynamic_cast<const GCTAEventCube*>(events()) == NULL) {
1276  std::string msg = "Specified file \""+cntcube.url()+"\" is not a CTA "
1277  "counts cube. Please provide a counts cube.";
1279  }
1280 
1281  // Load exposure cube
1282  GCTACubeExposure exposure(expcube);
1283 
1284  // Load point spread function cube
1285  GCTACubePsf psf(psfcube);
1286 
1287  // Load energy dispersion cube
1288  GCTACubeEdisp edisp(edispcube);
1289 
1290  // Load background cube
1291  GCTACubeBackground background(bkgcube);
1292 
1293  // Set exposure cube, point spread function cube, energy dispersion
1294  // cube, and background cube as response for this observation
1295  response(exposure, psf, edisp, background);
1296 
1297  // Return
1298  return;
1299 }
1300 
1301 
1302 /***********************************************************************//**
1303  * @brief Save CTA observation into FITS file.
1304  *
1305  * @param[in] filename FITS filename.
1306  * @param[in] clobber Overwrite existing FITS file? (default: false)
1307  *
1308  * Saves the CTA observation into a FITS file.
1309  *
1310  * If the CTA observation contains an event list, the method will write an
1311  * events and a Good Time Intervals extension to the FITS file. The names
1312  * for both extension can be optionally specified in the filename using
1313  * the format
1314  *
1315  * filename[EVENTS;GTI]
1316  *
1317  * where the string in the squared bracket defines the extension names. The
1318  * part before the semi-colon is the events extension name and the part after
1319  * the semi-colon is the Good Time Intervals extension name.
1320  *
1321  * If the CTA observation contains an event cube, the method will write the
1322  * cube into the primary image, followed by binary tables containing the
1323  * energy boundaries and the Good Time Intervals. The extension names of
1324  * these binary tables are `EBOUNDS` and `GTI`, and cannot be modified.
1325  ***************************************************************************/
1326 void GCTAObservation::save(const GFilename& filename,
1327  const bool& clobber) const
1328 {
1329  // Set default events and Good Time Intervals extension name and extract
1330  // a possible overwrite from the extension name argument of the filename.
1331  // The specific format that is implemented is [events;gti], where the
1332  // part before the semi-colon is the events extension name and the part
1333  // after the semi-colon is the Good Time Intervals extension name.
1334  std::string evtname = gammalib::extname_cta_events;
1335  std::string gtiname = gammalib::extname_gti;
1336  if (filename.has_extname()) {
1337  std::vector<std::string> extnames = gammalib::split(filename.extname(), ";");
1338  if (extnames.size() > 0) {
1339  evtname = gammalib::strip_whitespace(extnames[0]);
1340  }
1341  if (extnames.size() > 1) {
1342  gtiname = gammalib::strip_whitespace(extnames[1]);
1343  }
1344  }
1345 
1346  // Open or create FITS file. Since we accept here a special structure
1347  // for the extension that cfitsio does not understand we only pass
1348  // the URL without extension name to the GFits constructor.
1349  GFits fits(filename.url(), true);
1350 
1351  // Write data into FITS file
1352  write(fits, evtname, gtiname);
1353 
1354  // Save FITS file
1355  fits.save(clobber);
1356 
1357  // Close FITS file
1358  fits.close();
1359 
1360  // Return
1361  return;
1362 }
1363 
1364 
1365 /***********************************************************************//**
1366  * @brief Return event type string
1367  *
1368  * @return Event type string.
1369  *
1370  * Returns "EventList" if the observation contains an event list,
1371  * "CountsCube" if it contains a counts cube, "Events" if it container an
1372  * unknown type of events (which should never occur), and an empty string if
1373  * no events have been allocated.
1374  ***************************************************************************/
1375 std::string GCTAObservation::eventtype(void) const
1376 {
1377  // Initialise empty event type string
1378  std::string eventtype = "";
1379 
1380  // Continue only if events are allocated
1381  if (m_events != NULL) {
1382 
1383  // Case A: we have a list
1384  GCTAEventList* list = dynamic_cast<GCTAEventList*>(m_events);
1385  if (list != NULL) {
1386  eventtype = "EventList";
1387  }
1388 
1389  // Case B: we have a cube
1390  else {
1391  GCTAEventCube* cube = dynamic_cast<GCTAEventCube*>(m_events);
1392  if (cube != NULL) {
1393  eventtype = "CountsCube";
1394  }
1395 
1396  // Case C: we don't know what we have
1397  else {
1398  eventtype = "Events";
1399  }
1400  }
1401 
1402  } // endif: events were allocated
1403 
1404  // Return event type
1405  return eventtype;
1406 }
1407 
1408 
1409 /***********************************************************************//**
1410  * @brief Dispose events
1411  *
1412  * Disposes the events to save memory. The method only applies to event
1413  * lists. If does nothing in case that the observation does not contain an
1414  * event list.
1415  ***************************************************************************/
1417 {
1418  // Get pointer to event list
1419  GCTAEventList* list = dynamic_cast<GCTAEventList*>(m_events);
1420 
1421  // Dispose event list if pointer is valid
1422  if (list != NULL) {
1423  list->dispose();
1424  }
1425 
1426  // Return
1427  return;
1428 }
1429 
1430 
1431 /***********************************************************************//**
1432  * @brief Print CTA observation information
1433  *
1434  * @param[in] chatter Chattiness.
1435  * @return String containing observation information.
1436  ***************************************************************************/
1437 std::string GCTAObservation::print(const GChatter& chatter) const
1438 {
1439  // Initialise result string
1440  std::string result;
1441 
1442  // Continue only if chatter is not silent
1443  if (chatter != SILENT) {
1444 
1445  // Append header
1446  result.append("=== GCTAObservation ===");
1447 
1448  // Append information
1449  result.append("\n"+gammalib::parformat("Name")+name());
1450  result.append("\n"+gammalib::parformat("Identifier")+id());
1451  result.append("\n"+gammalib::parformat("Instrument")+instrument());
1452  result.append("\n"+gammalib::parformat("Event file")+eventfile().url());
1453  result.append("\n"+gammalib::parformat("Event type")+eventtype());
1454  result.append("\n"+gammalib::parformat("Statistic")+statistic());
1455  result.append("\n"+gammalib::parformat("Ontime"));
1456  result.append(gammalib::str(ontime())+" s");
1457  result.append("\n"+gammalib::parformat("Livetime"));
1458  result.append(gammalib::str(livetime())+" s");
1459  result.append("\n"+gammalib::parformat("Deadtime correction"));
1460  result.append(gammalib::str(m_deadc));
1461 
1462  // Append user energy threshold information
1463  result.append("\n"+gammalib::parformat("User energy range"));
1464  if (m_lo_user_thres > 0.0 && m_hi_user_thres) {
1465  result.append(gammalib::str(m_lo_user_thres));
1466  result.append(" - ");
1467  result.append(gammalib::str(m_hi_user_thres));
1468  result.append(" TeV");
1469  }
1470  else if (m_lo_user_thres > 0.0) {
1471  result.append("> ");
1472  result.append(gammalib::str(m_lo_user_thres));
1473  result.append(" TeV");
1474  }
1475  else if (m_hi_user_thres > 0.0) {
1476  result.append("< ");
1477  result.append(gammalib::str(m_hi_user_thres));
1478  result.append(" TeV");
1479  }
1480  else {
1481  result.append("undefined");
1482  }
1483 
1484  // Append detailed information
1485  GChatter reduced_chatter = gammalib::reduce(chatter);
1486  if (reduced_chatter > SILENT) {
1487 
1488  // Append pointing
1489  result.append("\n"+pointing().print(reduced_chatter));
1490 
1491  // Append response
1492  if (m_response != NULL) {
1493  result.append("\n"+m_response->print(reduced_chatter));
1494  }
1495  else {
1496  result.append("\n"+gammalib::parformat("Response function"));
1497  result.append("undefined");
1498  }
1499 
1500  // Append events
1501  if (m_events != NULL) {
1502  result.append("\n"+m_events->print(reduced_chatter));
1503  }
1504 
1505  // Append Off regions
1506  result.append("\n"+off_regions().print(reduced_chatter));
1507 
1508  } // endif: appended detailed information
1509 
1510  } // endif: chatter was not silent
1511 
1512  // Return result
1513  return result;
1514 }
1515 
1516 
1517 /*==========================================================================
1518  = =
1519  = Private methods =
1520  = =
1521  ==========================================================================*/
1522 
1523 /***********************************************************************//**
1524  * @brief Initialise class members
1525  ***************************************************************************/
1527 {
1528  // Initialise members
1529  m_instrument = "CTA";
1530  m_object.clear();
1531  m_eventfile.clear();
1532  m_bgdfile.clear();
1533  m_response = NULL;
1534  m_pointing.clear();
1535  m_off_regions.clear();
1536  m_ontime = 0.0;
1537  m_livetime = 0.0;
1538  m_deadc = 1.0;
1539  m_ra_obj = 0.0;
1540  m_dec_obj = 0.0;
1541  m_lo_user_thres = 0.0;
1542  m_hi_user_thres = 0.0;
1543  m_n_tels = 0;
1544 
1545  // Return
1546  return;
1547 }
1548 
1549 
1550 /***********************************************************************//**
1551  * @brief Copy class members
1552  *
1553  * @param[in] obs CTA observation.
1554  ***************************************************************************/
1556 {
1557  // Copy members
1558  m_instrument = obs.m_instrument;
1559  m_object = obs.m_object;
1560  m_eventfile = obs.m_eventfile;
1561  m_bgdfile = obs.m_bgdfile;
1562  m_pointing = obs.m_pointing;
1564  m_ontime = obs.m_ontime;
1565  m_livetime = obs.m_livetime;
1566  m_deadc = obs.m_deadc;
1567  m_ra_obj = obs.m_ra_obj;
1568  m_dec_obj = obs.m_dec_obj;
1571  m_n_tels = obs.m_n_tels;
1572 
1573  // Clone members
1574  m_response = (obs.m_response != NULL) ? obs.m_response->clone() : NULL;
1575 
1576  // Return
1577  return;
1578 }
1579 
1580 
1581 /***********************************************************************//**
1582  * @brief Delete class members
1583  ***************************************************************************/
1585 {
1586  // Free memory
1587  if (m_response != NULL) delete m_response;
1588 
1589  // Initialise pointers
1590  m_response = NULL;
1591 
1592  // Return
1593  return;
1594 }
1595 
1596 
1597 /***********************************************************************//**
1598  * @brief Read observation attributes
1599  *
1600  * @param[in] hdu FITS HDU.
1601  *
1602  * Reads CTA observation attributes from HDU. Mandatory attributes are
1603  *
1604  * RA_PNT - Right Ascension of pointing
1605  * DEC_PNT - Declination of pointing
1606  * ONTIME - Exposure time
1607  * LIVETIME - Livetime
1608  *
1609  * and optional attributes are
1610  *
1611  * OBJECT - Name of observed object
1612  * DEADC - Deadtime correction
1613  * RA_OBJ - Right Ascension of observed object,
1614  * DEC_OBJ - Declination of observed object,
1615  * OBS_ID - Observation identifier
1616  * ALT_PNT - Altitude of pointing above horizon
1617  * AZ_PNT - Azimuth of pointing
1618  * TELESCOP - Telescope name
1619  * N_TELS - Number of telescopes
1620  *
1621  * Based on RA_PNT and DEC_PNT, the CTA pointing direction is set. Note that
1622  * DEADC is computed using DEADC=LIVETIME/ONTIME
1623  *
1624  * @todo The actual reader is a minimal reader to accomodate as many
1625  * different datasets as possible. Once the CTA data format is fixed
1626  * the reader should have more mandatory attributes.
1627  ***************************************************************************/
1629 {
1630  // Read mandatory attributes
1631  double ra_pnt = hdu.real("RA_PNT");
1632  double dec_pnt = hdu.real("DEC_PNT");
1633  m_ontime = (hdu.has_card("ONTIME")) ? hdu.real("ONTIME") : 0.0;
1634  m_livetime = (hdu.has_card("LIVETIME")) ? hdu.real("LIVETIME") : 0.0;
1635 
1636  // Read optional attributes
1637  m_object = (hdu.has_card("OBJECT")) ? hdu.string("OBJECT") : "unknown";
1638  m_deadc = (hdu.has_card("DEADC")) ? hdu.real("DEADC") : 0.0;
1639  m_ra_obj = (hdu.has_card("RA_OBJ")) ? hdu.real("RA_OBJ") : 0.0;
1640  m_dec_obj = (hdu.has_card("DEC_OBJ")) ? hdu.real("DEC_OBJ") : 0.0;
1641  std::string id = (hdu.has_card("OBS_ID")) ? hdu.string("OBS_ID") : "";
1642  double alt = (hdu.has_card("ALT_PNT")) ? hdu.real("ALT_PNT") : 0.0;
1643  double az = (hdu.has_card("AZ_PNT")) ? hdu.real("AZ_PNT") : 0.0;
1644  m_instrument = (hdu.has_card("TELESCOP")) ? hdu.string("TELESCOP") : "CTA";
1645  m_n_tels = (hdu.has_card("N_TELS")) ? hdu.integer("N_TELS") : 0;
1646 
1647  // Kluge: deal with "H.E.S.S." telescope name
1648  if (m_instrument == "H.E.S.S.") {
1649  m_instrument = "HESS";
1650  }
1651 
1652  // Set attributes
1653  this->id(id);
1654 
1655  // Kluge: compute DEADC from livetime and ontime instead of using the
1656  // keyword value as the original event lists had this values badly
1657  // assigned
1658  if (m_ontime > 0) {
1660  }
1661  else {
1662  m_deadc = 0.0;
1663  }
1664 
1665  // Set pointing information
1666  GSkyDir pnt;
1667  pnt.radec_deg(ra_pnt, dec_pnt);
1668  m_pointing.dir(pnt);
1669  m_pointing.zenith(90.0-alt);
1670  m_pointing.azimuth(az);
1671 
1672  // Return
1673  return;
1674 }
1675 
1676 
1677 /***********************************************************************//**
1678  * @brief Write observation attributes
1679  *
1680  * @param[in] hdu FITS HDU.
1681  ***************************************************************************/
1683 {
1684  // Get time reference
1685  GTimeReference timeref = events()->gti().reference();
1686 
1687  // Compute some attributes
1688  double ra_pnt = m_pointing.dir().ra_deg();
1689  double dec_pnt = m_pointing.dir().dec_deg();
1690  double alt = 90.0 - m_pointing.zenith();
1691  double az = m_pointing.azimuth();
1692  double tstart = events()->tstart().convert(timeref);
1693  double tstop = events()->tstop().convert(timeref);
1694  double telapse = events()->gti().telapse();
1695  double ontime = events()->gti().ontime();
1696  double deadc = (ontime > 0.0 && livetime() > 0.0) ?
1697  livetime() / ontime : m_deadc;
1698  std::string utc_obs = events()->tstart().utc();
1699  std::string utc_end = events()->tstop().utc();
1700  std::string date_obs = utc_obs.substr(0, 10);
1701  std::string time_obs = utc_obs.substr(11, 8);
1702  std::string date_end = utc_end.substr(0, 10);
1703  std::string time_end = utc_end.substr(11, 8);
1704 
1705  // Set observation information
1706  hdu.card("CREATOR", "GammaLib", "Program which created the file");
1707  hdu.card("TELESCOP", instrument(), "Telescope");
1708  hdu.card("OBS_ID", id(), "Observation identifier");
1709  hdu.card("DATE-OBS", date_obs, "Observation start date");
1710  hdu.card("TIME-OBS", time_obs, "Observation start time");
1711  hdu.card("DATE-END", date_end, "Observation end date");
1712  hdu.card("TIME-END", time_end, "Observation end time");
1713 
1714  // Set observation time information
1715  hdu.card("TSTART", tstart, "[s] Mission time of start of observation");
1716  hdu.card("TSTOP", tstop, "[s] Mission time of end of observation");
1717  timeref.write(hdu);
1718  hdu.card("TELAPSE", telapse, "[s] Mission elapsed time");
1719  hdu.card("ONTIME", ontime, "[s] Total good time including deadtime");
1720  hdu.card("LIVETIME", livetime(), "[s] Total livetime");
1721  hdu.card("DEADC", deadc, "Deadtime correction factor");
1722  hdu.card("TIMEDEL", 1.0, "Time resolution");
1723 
1724  // Set pointing information
1725  hdu.card("OBJECT", object(), "Observed object");
1726  hdu.card("RA_OBJ", ra_obj(), "[deg] Target Right Ascension");
1727  hdu.card("DEC_OBJ", dec_obj(), "[deg] Target Declination");
1728  hdu.card("RA_PNT", ra_pnt, "[deg] Pointing Right Ascension");
1729  hdu.card("DEC_PNT", dec_pnt, "[deg] Pointing Declination");
1730  hdu.card("ALT_PNT", alt, "[deg] Average altitude of pointing");
1731  hdu.card("AZ_PNT", az, "[deg] Average azimuth of pointing");
1732  hdu.card("RADECSYS", "FK5", "Coordinate system");
1733  hdu.card("EQUINOX", 2000.0, "Epoch");
1734  hdu.card("CONV_DEP", 0.0, "Convergence depth of telescopes");
1735  hdu.card("CONV_RA", 0.0, "[deg] Convergence Right Ascension");
1736  hdu.card("CONV_DEC", 0.0, "[deg] Convergence Declination");
1737  hdu.card("OBSERVER", "string", "Observer");
1738 
1739  // Telescope information
1740  hdu.card("N_TELS", n_tels(), "Number of telescopes in event list");
1741  hdu.card("TELLIST", "string", "Telescope IDs");
1742  hdu.card("GEOLAT", 0.0, "[deg] Geographic latitude of array centre");
1743  hdu.card("GEOLON", 0.0, "[deg] Geographic longitude of array centre");
1744  hdu.card("ALTITUDE", 0.0, "[km] Altitude of array centre");
1745 
1746  // Other information
1747  hdu.card("EUNIT", "TeV", "Energy unit");
1748  hdu.card("EVTVER", "draft1", "Event list version number");
1749 
1750  // Return
1751  return;
1752 }
void clear(void)
Clear CTA pointing.
const std::string & statistic(void) const
Return optimizer statistic.
const GSkyRegions & off_regions(void) const
Return sky off regions.
GFitsHDU * at(const int &extno)
Get pointer to HDU.
Definition: GFits.cpp:233
#define G_GTI
virtual ~GCTAObservation(void)
Destructor.
CTA exposure cube class.
const GFilename & eventfile(void) const
Return event file name.
void clear(void)
Clear object.
CTA event list class interface definition.
bool contains(const int &extno) const
Check if HDU exists in FITS file.
Definition: GFits.hpp:282
bool has_card(const int &cardno) const
Check existence of header card.
Definition: GFitsHDU.hpp:233
double dec_deg(void) const
Returns Declination in degrees.
Definition: GSkyDir.hpp:256
void xml_check_par(const std::string &origin, const std::string &name, const int &number)
Checks whether a parameter has occured once.
Definition: GTools.cpp:1855
void save(const GFilename &filename, const bool &clobber=false) const
Save CTA observation into FITS file.
CTA cube-style response function class definition.
double m_ontime
Ontime (seconds)
virtual void write(GFits &fits) const
Write CTA events and Good Time Intervals into FITS file.
virtual void roi(const GRoi &roi)
Set ROI.
void write(GFitsHDU &hdu) const
Write time reference into FITS header.
GFilename m_eventfile
Event filename.
void ebounds(const GEbounds &ebounds)
Set energy boundaries.
Definition: GEvents.cpp:136
Abstract FITS extension base class.
Definition: GFitsHDU.hpp:51
#define G_RESPONSE_GET
#define G_RESPONSE_SET
virtual void read(const GXmlElement &xml)
Read observation from XML element.
CTA cube background class.
GEvents * m_events
Pointer to event container.
GEbounds ebounds(void) const
Get energy boundaries.
const GFilename & filename(void) const
Return FITS filename.
Definition: GFits.hpp:313
void load(const GFilename &filename)
Load unbinned or binned analysis data.
virtual std::string print(const GChatter &chatter=NORMAL) const
Print CTA observation information.
CTA cube-style response function class.
const GCTAPointing & pointing(void) const
Return CTA pointing.
bool is_empty(void) const
Signal if filename is empty.
Definition: GFilename.hpp:160
CTA instrument response function class definition.
const GTime & tstop(void) const
Return stop time.
Definition: GEvents.hpp:158
CTA event list class.
XML element node class.
Definition: GXmlElement.hpp:48
void gti(const GGti &gti)
Set Good Time Intervals.
Definition: GEvents.cpp:154
bool has_extname(void) const
Signal if filename has an extension name.
Definition: GFilename.hpp:255
const GCTAObservation g_obs_fact_seed(true,"FACT")
const double & zenith(void) const
Return pointing zenith angle.
void dispose_events(void)
Dispose events.
std::string extname(const std::string &defaultname="") const
Return extension name.
Definition: GFilename.cpp:385
const GCTAObservation g_obs_astri_seed(true,"ASTRI")
std::vector< std::string > split(const std::string &s, const std::string &sep)
Split string.
Definition: GTools.cpp:983
const int & n_tels(void) const
Return number of telescopes.
Interface for the CTA region of interest class.
Definition: GCTARoi.hpp:49
void read(const GFitsTable &table)
Read energy boundaries from FITS table.
Definition: GEbounds.cpp:761
virtual int elements(void) const
Return number of GXMLElement children of node.
Definition: GXmlNode.cpp:586
Gammalib tools definition.
void read(const GXmlElement &xml)
Read pointing from XML element.
FITS file class.
Definition: GFits.hpp:63
const std::string extname_gti
Definition: GGti.hpp:44
FITS file class interface definition.
double m_lo_user_thres
User defined lower energy threshold.
void write(GXmlElement &xml) const
Write pointing information into XML element.
virtual GCTAResponse * clone(void) const =0
Clones object.
virtual double deadc(const GTime &time=GTime()) const
Return deadtime correction factor.
double m_hi_user_thres
User defined upper energy boundary.
std::string strip_whitespace(const std::string &arg)
Strip leading and trailing whitespace from string.
Definition: GTools.cpp:80
GFilename xml_file_reduce(const GXmlElement &xml, const std::string &filename)
Reduce file name provided for writing as XML attribute.
Definition: GTools.cpp:1946
Good time interval class interface definition.
void free_members(void)
Delete class members.
virtual double ontime(void) const
Return ontime.
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
#define G_LOAD2
CTA event bin container class.
virtual double livetime(void) const
Return livetime.
Calibration database class.
Definition: GCaldb.hpp:66
const GFilename & filename(void) const
Return regions file name.
const std::string & id(void) const
Return observation identifier.
Energy boundaries container class.
Definition: GEbounds.hpp:60
void remove(const int &extno)
Remove HDU from FITS file.
Definition: GFits.cpp:862
const GXmlAttribute * attribute(const int &index) const
Return attribute.
GCTARoi roi(void) const
Get Region of Interest.
const GCTAResponse * cta_rsp(const std::string &origin, const GResponse &rsp)
Retrieve CTA response from generic observation.
virtual void clear(void)
Clear CTA observation.
double m_deadc
Deadtime correction (livetime/ontime)
Filename class.
Definition: GFilename.hpp:62
#define G_ROI
GGti gti(void) const
Get Good Time Intervals.
const double & dec_obj(void) const
Return CTA object Declination.
std::string m_object
Object name.
const double & ra_obj(void) const
Return CTA object Right Ascension.
const std::string extname_cta_events
GXmlElement * xml_need_par(const std::string &origin, GXmlElement &xml, const std::string &name)
Return pointer to parameter with given name in XML element.
Definition: GTools.cpp:1637
virtual std::string instrument(void) const
Return instrument name.
Interface definition for the observation registry class.
const std::string & name(void) const
Return observation name.
Observation registry class definition.
void dispose(void) const
Dispose events.
double ra_deg(void) const
Returns Right Ascension in degrees.
Definition: GSkyDir.hpp:229
const GCTAObservation g_obs_magic_seed(true,"MAGIC")
#define G_LOAD1
CTA energy dispersion for stacked analysis.
GChatter
Definition: GTypemaps.hpp:33
GCTAResponse * m_response
Pointer to instrument response functions.
GCTAObservation & operator=(const GCTAObservation &obs)
Assignment operator.
double deadc(void) const
Return deadtime correction.
int integer(const std::string &keyname) const
Return card value as integer.
Definition: GFitsHDU.hpp:436
Good Time Interval class.
Definition: GGti.hpp:62
virtual void write(GXmlElement &xml) const
Write observation into XML element.
void init_members(void)
Initialise class members.
Abstract observation base class.
void radec_deg(const double &ra, const double &dec)
Set equatorial sky direction (degrees)
Definition: GSkyDir.cpp:224
void copy_members(const GCTAObservation &obs)
Copy class members.
virtual const GCTAResponse * response(void) const
Return pointer to CTA response function.
CTA region of interest class interface definition.
#define G_READ
CTA observation class interface definition.
virtual void read(const GXmlElement &xml)=0
void load(const std::string &rspname)
Load CTA response.
bool xml_has_par(const GXmlElement &xml, const std::string &name)
Checks if parameter with given name in XML element exists.
Definition: GTools.cpp:1596
CTA instrument response function class.
CTA instrument response function class.
Sky region container class.
Definition: GSkyRegions.hpp:56
const double & azimuth(void) const
Return pointing azimuth angle.
const GCTAObservation g_obs_veritas_seed(true,"VERITAS")
const GCTAObservation g_obs_hess_seed(true,"HESS")
std::string url(void) const
Return Uniform Resource Locator (URL)
Definition: GFilename.hpp:189
virtual GXmlElement * element(const int &index)
Return pointer to GXMLElement child.
Definition: GXmlNode.cpp:640
void free_members(void)
Delete class members.
void clear(void)
Clear file name.
Definition: GFilename.cpp:188
const GCTAObservation g_obs_cta_seed(true,"CTA")
CTA point spread function for cube analysis.
Definition: GCTACubePsf.hpp:62
void init_members(void)
Initialise class members.
void write_attributes(GFitsHDU &hdu) const
Write observation attributes.
void read(const GFitsTable &table)
Read Good Time Intervals and time reference from FITS table.
Definition: GGti.cpp:753
double m_dec_obj
Declination of object (degrees)
GChatter reduce(const GChatter &chatter)
Reduce chattiness by one level.
Definition: GTypemaps.hpp:65
virtual std::string print(const GChatter &chatter=NORMAL) const =0
Print content of object.
const double & ontime(void) const
Return ontime.
CTA event bin container class interface definition.
#define G_EBOUNDS
double m_ra_obj
Right Ascension of object (degrees)
const double & livetime(void) const
Return livetime.
virtual GObservation & operator=(const GObservation &obs)
Assignment operator.
virtual GEvents * events(void)
Return events.
std::string string(const std::string &keyname) const
Return card value as string.
Definition: GFitsHDU.hpp:410
std::string eventtype(void) const
Return event type string.
Exception handler interface definition.
std::string toupper(const std::string &s)
Convert string to upper case.
Definition: GTools.cpp:941
Implements a time reference.
GCTAPointing m_pointing
Pointing direction.
void read_attributes(const GFitsHDU &hdu)
Read observation attributes.
const GSkyDir & dir(void) const
Return pointing sky direction.
virtual GCTAObservation * clone(void) const
Clone CTA observation.
virtual void write(GXmlElement &xml) const =0
void read(const GXmlElement &xml)
Read region of interest from XML element.
Definition: GCTARoi.cpp:265
virtual void write(GFits &file) const
Write CTA event cube into FITS file.
CTA observation class.
Abstract instrument response base class.
Definition: GResponse.hpp:77
std::string parformat(const std::string &s, const int &indent=0)
Convert string in parameter format.
Definition: GTools.cpp:1143
const GTime & tstart(void) const
Return start time.
Definition: GEvents.hpp:146
Integration class interface definition.
GFitsHeaderCard & card(const int &cardno)
Return header card.
Definition: GFitsHDU.hpp:259
Sky direction class.
Definition: GSkyDir.hpp:62
#define G_WRITE
void close(void)
Close FITS file.
Definition: GFits.cpp:1342
GSkyRegions m_off_regions
Off regions.
std::string utc(const int &precision=0) const
Return time as string in UTC time system.
Definition: GTime.cpp:465
void caldb(const GCaldb &caldb)
Set calibration database.
const GXmlElement * xml_get_par(const std::string &origin, const GXmlElement &xml, const std::string &name)
Return pointer to parameter with given name in XML element.
Definition: GTools.cpp:1689
const double & ontime(void) const
Returns ontime.
Definition: GGti.hpp:240
virtual void read(const GFits &fits)
Read events from FITS file.
virtual void read(const GFits &file)
Read CTA event cube from FITS file.
double convert(const GTimeReference &ref) const
Return time in specified reference.
Definition: GTime.cpp:698
std::string m_instrument
Instrument name.
Filename class interface definition.
virtual std::string print(const GChatter &chatter=NORMAL) const =0
Print content of object.
int m_n_tels
Number of telescopes.
double todouble(const std::string &arg)
Convert string into double precision value.
Definition: GTools.cpp:926
GCTAObservation(void)
Void constructor.
double m_livetime
Livetime (seconds)
GFilename xml_file_expand(const GXmlElement &xml, const std::string &filename)
Expand file name provided as XML attribute for loading.
Definition: GTools.cpp:1889
std::string m_bgdfile
Background filename.
std::string str(const unsigned short int &value)
Convert unsigned short integer value into string.
Definition: GTools.cpp:489