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