GammaLib  2.0.0
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
GCTAEventCube.cpp
Go to the documentation of this file.
1 /***************************************************************************
2  * GCTAEventCube.cpp - CTA event bin container class *
3  * ----------------------------------------------------------------------- *
4  * copyright (C) 2010-2021 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 GCTAEventCube.cpp
23  * @brief CTA event bin container class implementation
24  * @author Juergen Knoedlseder
25  */
26 
27 /* __ Includes ___________________________________________________________ */
28 #ifdef HAVE_CONFIG_H
29 #include <config.h>
30 #endif
31 #include "GException.hpp"
32 #include "GTools.hpp"
33 #include "GFits.hpp"
34 #include "GCTAEventCube.hpp"
35 #include "GCTATypemaps.hpp"
36 
37 /* __ Method name definitions ____________________________________________ */
38 #define G_NAXIS "GCTAEventCube::naxis(int)"
39 #define G_ENERGY "GCTAEventCube::energy(int&)"
40 #define G_SET_DIRECTIONS "GCTAEventCube::set_directions()"
41 #define G_SET_ENERGIES "GCTAEventCube::set_energies()"
42 #define G_SET_TIMES "GCTAEventCube::set_times()"
43 #define G_SET_BIN "GCTAEventCube::set_bin(int&)"
44 
45 /* __ Macros _____________________________________________________________ */
46 
47 /* __ Coding definitions _________________________________________________ */
48 
49 /* __ Debug definitions __________________________________________________ */
50 
51 
52 
53 /*==========================================================================
54  = =
55  = Constructors/destructors =
56  = =
57  ==========================================================================*/
58 
59 /***********************************************************************//**
60  * @brief Void constructor
61  ***************************************************************************/
63 {
64  // Initialise members
65  init_members();
66 
67  // Return
68  return;
69 }
70 
71 
72 /***********************************************************************//**
73  * @brief File name constructor
74  *
75  * @param[in] filename Counts cube filename.
76  *
77  * Constructs instance of events cube from a counts cube file.
78  ***************************************************************************/
80 {
81  // Initialise members
82  init_members();
83 
84  // Load counts cube
85  load(filename);
86 
87  // Return
88  return;
89 }
90 
91 
92 /***********************************************************************//**
93  * @brief Constructor
94  *
95  * @param[in] map Sky map.
96  * @param[in] ebds Energy boundaries.
97  * @param[in] gti Good Time intervals.
98  *
99  * Constructs instance of events cube from a sky map, energy boundaries and
100  * Good Time Intervals. All event cube weights are set to unity.
101  ***************************************************************************/
103  const GEbounds& ebds,
104  const GGti& gti) : GEventCube()
105 {
106  // Initialise members
107  init_members();
108 
109  // Set sky map, energy boundaries and GTI
110  m_map = map;
111  this->ebounds(ebds);
112  this->gti(gti);
113 
114  // Set sky directions
115  set_directions();
116 
117  // Set energies
118  set_energies();
119 
120  // Set times
121  set_times();
122 
123  // Set all weights to unity
124  m_weights = m_map;
125  m_weights = 1.0;
126 
127  // Return
128  return;
129 }
130 
131 
132 /***********************************************************************//**
133  * @brief Constructor
134  *
135  * @param[in] map Sky map.
136  * @param[in] weights Event cube weights.
137  * @param[in] ebds Energy boundaries.
138  * @param[in] gti Good Time intervals.
139  *
140  * Constructs instance of events cube from a sky map, a map of weights,
141  * energy boundaries and Good Time Intervals.
142  ***************************************************************************/
144  const GSkyMap& weights,
145  const GEbounds& ebds,
146  const GGti& gti) : GEventCube()
147 {
148  // Initialise members
149  init_members();
150 
151  // Set sky map, energy boundaries and GTI
152  m_map = map;
153  this->ebounds(ebds);
154  this->gti(gti);
155 
156  // Set weight map
157  m_weights = weights;
158 
159  // Set sky directions
160  set_directions();
161 
162  // Set energies
163  set_energies();
164 
165  // Set times
166  set_times();
167 
168  // Return
169  return;
170 }
171 
172 
173 /***********************************************************************//**
174  * @brief Copy constructor
175  *
176  * @param[in] cube Event cube.
177  ***************************************************************************/
179 {
180  // Initialise members
181  init_members();
182 
183  // Copy members
184  copy_members(cube);
185 
186  // Return
187  return;
188 }
189 
190 
191 /***********************************************************************//**
192  * @brief Destructor
193  ***************************************************************************/
195 {
196  // Free members
197  free_members();
198 
199  // Return
200  return;
201 }
202 
203 
204 /*==========================================================================
205  = =
206  = Operators =
207  = =
208  ==========================================================================*/
209 
210 /***********************************************************************//**
211  * @brief Assignment operator
212  *
213  * @param[in] cube Event cube.
214  * @return Event cube
215  ***************************************************************************/
217 {
218  // Execute only if object is not identical
219  if (this != &cube) {
220 
221  // Copy base class members
222  this->GEventCube::operator=(cube);
223 
224  // Free members
225  free_members();
226 
227  // Initialise members
228  init_members();
229 
230  // Copy members
231  copy_members(cube);
232 
233  } // endif: object was not identical
234 
235  // Return this object
236  return *this;
237 }
238 
239 
240 /***********************************************************************//**
241  * @brief Event bin access operator
242  *
243  * @param[in] index Event index [0,...,size()-1].
244  *
245  * Returns pointer to an event bin.
246  ***************************************************************************/
248 {
249  // Set event bin
250  set_bin(index);
251 
252  // Return pointer
253  return (&m_bin);
254 }
255 
256 
257 /***********************************************************************//**
258  * @brief Event bin access operator (const version)
259  *
260  * @param[in] index Event index [0,...,size()-1].
261  *
262  * Returns pointer to an event bin.
263  ***************************************************************************/
264 const GCTAEventBin* GCTAEventCube::operator[](const int& index) const
265 {
266  // Set event bin (circumvent const correctness)
267  (const_cast<GCTAEventCube*>(this))->set_bin(index);
268 
269  // Return pointer
270  return (&m_bin);
271 }
272 
273 
274 /*==========================================================================
275  = =
276  = Public methods =
277  = =
278  ==========================================================================*/
279 
280 /***********************************************************************//**
281  * @brief Clear CTA event cube
282  ***************************************************************************/
284 {
285  // Free class members (base and derived classes, derived class first)
286  free_members();
287  this->GEventCube::free_members();
288  this->GEvents::free_members();
289 
290  // Initialise members
291  this->GEvents::init_members();
292  this->GEventCube::init_members();
293  init_members();
294 
295  // Return
296  return;
297 }
298 
299 
300 /***********************************************************************//**
301  * @brief Clone CTA event cube
302  *
303  * @return Pointer to deep copy of CTA event cube.
304  ***************************************************************************/
306 {
307  return new GCTAEventCube(*this);
308 }
309 
310 
311 /***********************************************************************//**
312  * @brief Return number of bins in event cube
313  ***************************************************************************/
314 int GCTAEventCube::size(void) const
315 {
316  // Compute number of bins
317  int nbins = m_map.npix() * m_map.nmaps();
318 
319  // Return number of bins
320  return nbins;
321 }
322 
323 
324 /***********************************************************************//**
325  * @brief Return number of bins in axis
326  *
327  * @param[in] axis Axis [0,...,dim()-1]
328  * @return Number of bins in specified axis
329  *
330  * @exception GException::out_of_range
331  * Axis is out of range.
332  *
333  * Returns the number of bins along a given event cube @p axis.
334  ***************************************************************************/
335 int GCTAEventCube::naxis(const int& axis) const
336 {
337  // Optionally check if the axis is valid
338  #if defined(G_RANGE_CHECK)
339  if (axis < 0 || axis >= dim()) {
340  throw GException::out_of_range(G_NAXIS, "CTA event cube axis", axis, dim());
341  }
342  #endif
343 
344  // Set result
345  int naxis = 0;
346  switch (axis) {
347  case 0:
348  naxis = m_map.nx();
349  break;
350  case 1:
351  naxis = m_map.ny();
352  break;
353  case 2:
354  naxis = m_map.nmaps();
355  break;
356  }
357 
358  // Return result
359  return naxis;
360 }
361 
362 
363 /***********************************************************************//**
364  * @brief Load CTA event cube from FITS file
365  *
366  * @param[in] filename FITS filename.
367  *
368  * Loads the event cube from a FITS file. See the read() method for more
369  * information about the structure of the FITS file.
370  *
371  * The method clears the object before loading, thus any events residing in
372  * the object before loading will be lost.
373  ***************************************************************************/
374 void GCTAEventCube::load(const GFilename& filename)
375 {
376  // Clear object
377  clear();
378 
379  // Open counts map FITS file
380  GFits fits(filename);
381 
382  // Load counts map
383  read(fits);
384 
385  // Close FITS file
386  fits.close();
387 
388  // Return
389  return;
390 }
391 
392 
393 /***********************************************************************//**
394  * @brief Save CTA event cube into FITS file
395  *
396  * @param[in] filename FITS filename.
397  * @param[in] clobber Overwrite existing FITS file.
398  *
399  * Save the CTA event cube into FITS file. See the write() method for more
400  * information about the structure of the FITS file.
401  ***************************************************************************/
402 void GCTAEventCube::save(const GFilename& filename,
403  const bool& clobber) const
404 {
405  // Create empty FITS file
406  GFits fits;
407 
408  // Write event cube
409  write(fits);
410 
411  // Save FITS file
412  fits.saveto(filename, clobber);
413 
414  // Return
415  return;
416 }
417 
418 
419 /***********************************************************************//**
420  * @brief Read CTA event cube from FITS file
421  *
422  * @param[in] fits FITS file.
423  *
424  * Read an event cube from a FITS file. The following HDUs will be read
425  *
426  * COUNTS - Counts cube (or primary extension if COUNTS does not exist)
427  * WEIGHTS - Weights for each counts cube bin (optional)
428  * EBOUNDS - Energy boundaries
429  * GTI - Good Time Intervals
430  *
431  * The method clears the event cube before reading, thus any events residing
432  * in the event cube will be lost.
433  ***************************************************************************/
434 void GCTAEventCube::read(const GFits& fits)
435 {
436  // Clear object
437  clear();
438 
439  // Set counts cube HDU
440  std::string counts_hdu("Primary");
442  counts_hdu = gammalib::extname_cta_counts;
443  }
444 
445  // Get HDUs
446  const GFitsImage& hdu_cntmap = *fits.image(counts_hdu);
447  const GFitsTable& hdu_ebounds = *fits.table(gammalib::extname_ebounds);
448  const GFitsTable& hdu_gti = *fits.table(gammalib::extname_gti);
449 
450  // Load counts map
451  read_cntmap(hdu_cntmap);
452 
453  // Load energy boundaries
454  read_ebds(hdu_ebounds);
455 
456  // Load GTIs
457  read_gti(hdu_gti);
458 
459  // If a WEIGHTS HDU exist then load it from the FITS file ...
461 
462  // Get WEIGHTS HDU
463  const GFitsImage& hdu_weights = *fits.image(gammalib::extname_cta_weights);
464 
465  // Read HDU
466  m_weights.read(hdu_weights);
467 
468  }
469 
470  // ... otherwise set the weight map to unity
471  else {
472  m_weights = m_map;
473  m_weights = 1.0;
474  }
475 
476  // Return
477  return;
478 }
479 
480 
481 /***********************************************************************//**
482  * @brief Write CTA event cube into FITS file.
483  *
484  * @param[in] fits FITS file.
485  *
486  * Writes CTA event cube into a FITS file. The following HDUs will be written
487  *
488  * COUNTS - Counts cube
489  * WEIGHTS - Weights for each counts cube bin
490  * EBOUNDS - Energy boundaries
491  * GTI - Good Time Intervals
492  *
493  * The counts cube will be written as a double precision sky map. The
494  * weighing cube contains the weight for each counts cube, computed as the
495  * fraction of the event bin that has been selected in filling the counts
496  * cube. This allows to account for proper stacking of observations with
497  * different energy thresholds, or different regions of interest. The energy
498  * boundaries for all counts cube layers are also written, as well as the
499  * Good Time Intervals that have been used in generating the counts cube.
500  ***************************************************************************/
501 void GCTAEventCube::write(GFits& fits) const
502 {
503  // Remove HDUs if they exist already
504  if (fits.contains(gammalib::extname_gti)) {
506  }
509  }
512  }
515  }
516  if (fits.contains("Primary")) {
517  fits.remove("Primary");
518  }
519 
520  // Write counts cube
522 
523  // Write cube weighting
525 
526  // Write energy boundaries
527  ebounds().write(fits);
528 
529  // Write Good Time intervals
530  gti().write(fits);
531 
532  // Return
533  return;
534 }
535 
536 
537 /***********************************************************************//**
538  * @brief Return number of events in cube
539  *
540  * @return Number of events in cube, rounded to nearest integer value.
541  *
542  * Returns the total number of events in the cube, rounded to the nearest
543  * integer value. All cube bins with a negative content will be excluded
544  * from the sum.
545  ***************************************************************************/
546 int GCTAEventCube::number(void) const
547 {
548  // Initialise result
549  double number = 0.0;
550 
551  // Get pointer on skymap pixels
552  const double* pixels = m_map.pixels();
553 
554  // Sum event cube
555  if (size() > 0 && pixels != NULL) {
556  for (int i = 0; i < size(); ++i) {
557  if (pixels[i] > 0.0) {
558  number += pixels[i];
559  }
560  }
561  }
562 
563  // Return
564  return int(number+0.5);
565 }
566 
567 
568 /***********************************************************************//**
569  * @brief Set event cube counts from sky map
570  *
571  * @param[in] counts Event cube counts sky map.
572  *
573  * Sets event cube counts from sky map. The methods also sets all weights to
574  * unity.
575  ***************************************************************************/
576 void GCTAEventCube::counts(const GSkyMap& counts)
577 {
578  // Store sky map
579  m_map = counts;
580 
581  // Compute sky directions
582  set_directions();
583 
584  // If event cube has pointing then set DETX and DETY coordinates of
585  // instrument direction
586  if (m_has_pnt) {
587  set_detxy(m_pnt);
588  }
589 
590  // Set all weights to unity
591  m_weights = m_map;
592  m_weights = 1.0;
593 
594  // Return
595  return;
596 }
597 
598 
599 /***********************************************************************//**
600  * @brief Return energy of cube layer
601  *
602  * @param[in] index Event cube layer index [0,...,ebins()-1]
603  * @return Energy of event cube layer
604  *
605  * @exception GException::out_of_range
606  * Layer index is out of range.
607  *
608  * Returns the energy of the event cube layer specified by @p index.
609  ***************************************************************************/
610 const GEnergy& GCTAEventCube::energy(const int& index) const
611 {
612  // Optionally check if the index is valid
613  #if defined(G_RANGE_CHECK)
614  if (index < 0 || index >= ebins()) {
616  "CTA event cube energy axis", index, ebins());
617  }
618  #endif
619 
620  // Return energy
621  return (m_energies[index]);
622 }
623 
624 
625 /***********************************************************************//**
626  * @brief Print event cube information
627  *
628  * @param[in] chatter Chattiness.
629  * @return String containing event cube information.
630  ***************************************************************************/
631 std::string GCTAEventCube::print(const GChatter& chatter) const
632 {
633  // Initialise result string
634  std::string result;
635 
636  // Continue only if chatter is not silent
637  if (chatter != SILENT) {
638 
639  // Append header
640  result.append("=== GCTAEventCube ===");
641  result.append("\n"+gammalib::parformat("Number of events") +
642  gammalib::str(number()));
643  result.append("\n"+gammalib::parformat("Number of elements") +
644  gammalib::str(size()));
645  result.append("\n"+gammalib::parformat("Number of pixels") +
646  gammalib::str(npix()));
647  result.append("\n"+gammalib::parformat("Number of energy bins") +
648  gammalib::str(ebins()));
649 
650  // Append GTI interval
651  result.append("\n"+gammalib::parformat("Time interval"));
652  if (gti().size() > 0) {
653  result.append(gammalib::str(tstart().mjd()));
654  result.append(" - ");
655  result.append(gammalib::str(tstop().mjd())+" days");
656  }
657  else {
658  result.append("not defined");
659  }
660 
661  // Append pointing
662  result.append("\n"+gammalib::parformat("Pointing"));
663  if (m_has_pnt) {
664  result.append(m_pnt.dir().print());
665  }
666  else {
667  result.append("not defined");
668  }
669 
670  // Append energy intervals
671  if (gammalib::reduce(chatter) > SILENT) {
672  if (ebounds().size() > 0) {
673  result.append("\n"+ebounds().print(gammalib::reduce(chatter)));
674  }
675  else {
676  result.append("\n"+gammalib::parformat("Energy intervals") +
677  "not defined");
678  }
679  }
680  else {
681  result.append("\n"+gammalib::parformat("Energy interval"));
682  if (ebounds().size() > 0) {
683  result.append(gammalib::str(emin().TeV()));
684  result.append(" - ");
685  result.append(gammalib::str(emax().TeV())+" TeV");
686  }
687  else {
688  result.append("not defined");
689  }
690  }
691 
692  // Append skymap definition
693  if (gammalib::reduce(chatter) > SILENT) {
694  result.append("\n"+m_map.print(gammalib::reduce(chatter)));
695  }
696 
697  } // endif: chatter was not silent
698 
699  // Return result
700  return result;
701 }
702 
703 
704 /*==========================================================================
705  = =
706  = Private methods =
707  = =
708  ==========================================================================*/
709 
710 /***********************************************************************//**
711  * @brief Initialise class members
712  ***************************************************************************/
714 {
715  // Initialise members
716  m_map.clear();
717  m_weights.clear();
718  m_bin.clear();
719  m_time.clear();
720  m_pnt.clear();
721  m_has_pnt = false;
722  m_dirs.clear();
723  m_solidangle.clear();
724  m_energies.clear();
725  m_ewidth.clear();
726  m_ontime = 0.0;
727 
728  // Prepare event bin
729  init_bin();
730 
731  // Set CTA time reference for GTIs
732  m_gti.reference(GTimeReference(G_CTA_MJDREF, "s", "TT", "LOCAL"));
733 
734  // Return
735  return;
736 }
737 
738 
739 /***********************************************************************//**
740  * @brief Copy class members
741  *
742  * @param[in] cube Event cube.
743  *
744  * This method copies the class members from another event cube in the actual
745  * object. It also prepares the event bin member that will be returned in
746  * case of an operator access to the class.
747  ***************************************************************************/
749 {
750  // Copy members. Note that the event bin is not copied as it will
751  // be initialised later. The event bin serves just as a container of
752  // pointers, hence we do not want to copy over the pointers from the
753  // original class.
754  m_map = cube.m_map;
755  m_weights = cube.m_weights;
756  m_time = cube.m_time;
757  m_pnt = cube.m_pnt;
758  m_has_pnt = cube.m_has_pnt;
759  m_dirs = cube.m_dirs;
760  m_solidangle = cube.m_solidangle;
761  m_energies = cube.m_energies;
762  m_ewidth = cube.m_ewidth;
763  m_ontime = cube.m_ontime;
764 
765  // Copy GTIs
766  m_gti = cube.m_gti;
767 
768  // Prepare event bin
769  init_bin();
770 
771  // Return
772  return;
773 }
774 
775 
776 /***********************************************************************//**
777  * @brief Delete class members
778  ***************************************************************************/
780 {
781  // Return
782  return;
783 }
784 
785 
786 /***********************************************************************//**
787  * @brief Read CTA counts map from HDU.
788  *
789  * @param[in] hdu Image HDU.
790  *
791  * This method reads a CTA counts map from a FITS HDU. The counts map is
792  * stored in a GSkyMap object.
793  ***************************************************************************/
795 {
796  // Load counts map as sky map
797  m_map.read(hdu);
798 
799  // Set sky directions
800  set_directions();
801 
802  // If header contains pointing direction then also set DETX and DETY
803  // coordinates
804  if (hdu.has_card("RA_PNT") && hdu.has_card("DEC_PNT")) {
805 
806  // Read pointing direction
807  double ra_pnt = hdu.real("RA_PNT");
808  double dec_pnt = hdu.real("DEC_PNT");
809  GSkyDir pnt;
810  pnt.radec_deg(ra_pnt, dec_pnt);
811  m_pnt.dir(pnt);
812  m_has_pnt = true;
813 
814  // Set DETX and DETY coordinates of instrument direction
815  set_detxy(m_pnt);
816 
817  } // endif: header contained pointing direction
818 
819  // Return
820  return;
821 }
822 
823 
824 /***********************************************************************//**
825  * @brief Read energy boundaries from HDU.
826  *
827  * @param[in] hdu Energy boundaries table.
828  *
829  * Read the energy boundaries from the HDU.
830  ***************************************************************************/
832 {
833  // Read energy boundaries
834  m_ebounds.read(hdu);
835 
836  // Set log mean energies and energy widths
837  set_energies();
838 
839  // Return
840  return;
841 }
842 
843 
844 /***********************************************************************//**
845  * @brief Read GTIs from HDU.
846  *
847  * @param[in] hdu GTI table.
848  *
849  * Reads the Good Time Intervals from the HDU.
850  ***************************************************************************/
852 {
853  // Read Good Time Intervals
854  m_gti.read(hdu);
855 
856  // Set times
857  set_times();
858 
859  // Return
860  return;
861 }
862 
863 
864 /***********************************************************************//**
865  * @brief Set sky directions and solid angles of events cube.
866  *
867  * @exception GException::invalid_value
868  * No sky pixels found in event cube.
869  *
870  * This method computes the sky directions and solid angles for all event
871  * cube pixels. Sky directions are stored in an array of GCTAInstDir objects
872  * while solid angles are stored in units of sr in an array of double
873  * precision variables.
874  *
875  * A kluge has been introduced that handles invalid pixels. Invalid pixels
876  * may occur if a Hammer-Aitoff projection is used. In this case, pixels may
877  * lie outside the valid sky region. As invalid pixels lead to exceptions
878  * in the WCS classes, we simply need to catch the exceptions here. Invalid
879  * pixels are signaled by setting the solid angle of the pixel to 0.
880  ***************************************************************************/
882 {
883  // Throw an error if we have no sky pixels
884  if (npix() < 1) {
885  std::string msg = "CTA event cube contains no sky pixels. Please "
886  "provide a valid CTA event cube.";
888  }
889 
890  // Clear old pixel directions and solid angle
891  m_dirs.clear();
892  m_solidangle.clear();
893 
894  // Reserve space for pixel directions and solid angles
895  m_dirs.reserve(npix());
896  m_solidangle.reserve(npix());
897 
898  // Set pixel directions and solid angles
899  for (int iy = 0; iy < ny(); ++iy) {
900  for (int ix = 0; ix < nx(); ++ix) {
901  try {
902  GSkyPixel pixel = GSkyPixel(double(ix), double(iy));
903  m_dirs.push_back(GCTAInstDir(m_map.pix2dir(pixel)));
904  m_solidangle.push_back(m_map.solidangle(pixel));
905  }
906  catch (GException::invalid_argument& e) {
907  m_dirs.push_back(GCTAInstDir());
908  m_solidangle.push_back(0.0);
909  }
910  }
911  }
912 
913  // Return
914  return;
915 }
916 
917 
918 /***********************************************************************//**
919  * @brief Set DETX and DETY coordinates
920  *
921  * @param[in] pnt CTA pointing.
922  *
923  * Computes DETX and DETY coordinates for a given pointing direction.
924  ***************************************************************************/
926 {
927  // Get number of instrument directions
928  int size = m_dirs.size();
929 
930  // Loop over all intstrument directions
931  for (int i = 0; i < size; ++i) {
932 
933  // Get sky direction
934  GSkyDir skydir = m_dirs[i].dir();
935 
936  // Set instrument direction from sky direction
937  m_dirs[i] = pnt.instdir(skydir);
938 
939  } // endfor: looped over all instrument directions
940 
941  // Return
942  return;
943 }
944 
945 
946 /***********************************************************************//**
947  * @brief Set log mean energies and energy widths of event cube.
948  *
949  * @exception GException::invalid_value
950  * No energy boundaries found in event cube.
951  *
952  * This method computes the log mean energies and the energy widths of the
953  * event cube. The log mean energies and energy widths are stored unit
954  * independent in arrays of GEnergy objects.
955  ***************************************************************************/
957 {
958  // Determine number of energy bins
959  int ebins = ebounds().size();
960 
961  // Throw an error if we have no energy bins
962  if (ebins < 1) {
963  std::string msg = "CTA event cube contains no energy boundaries. "
964  "Please provide a valid CTA event cube.";
966  }
967 
968  // Clear old bin energies and energy widths
969  m_energies.clear();
970  m_ewidth.clear();
971 
972  // Reserve space for bin energies and energy widths
973  m_energies.reserve(ebins);
974  m_ewidth.reserve(ebins);
975 
976  // Setup bin energies and energy widths
977  for (int i = 0; i < ebins; ++i) {
978  m_energies.push_back(ebounds().elogmean(i));
979  m_ewidth.push_back(ebounds().emax(i) - ebounds().emin(i));
980  }
981 
982  // Return
983  return;
984 }
985 
986 
987 /***********************************************************************//**
988  * @brief Set mean event time and ontime of event cube
989  *
990  * @exception GException::invalid_value
991  * No Good Time Intervals found.
992  *
993  * Computes the mean time of the event cube by taking the mean between start
994  * and stop time. Computes also the ontime by summing up of all good time
995  * intervals.
996  *
997  * @todo Could add a more sophisticated mean event time computation that
998  * weights by the length of the GTIs, yet so far we do not really use
999  * the mean event time, hence there is no rush to implement this.
1000  ***************************************************************************/
1002 {
1003  // Throw an error if GTI is empty
1004  if (m_gti.size() < 1) {
1005  std::string msg = "No Good Time Intervals have been found in event "
1006  "cube. Every CTA event cube needs a definition "
1007  "of the Good Time Intervals.";
1009  }
1010 
1011  // Compute mean time
1012  m_time = m_gti.tstart() + 0.5 * (m_gti.tstop() - m_gti.tstart());
1013 
1014  // Set ontime
1015  m_ontime = m_gti.ontime();
1016 
1017  // Return
1018  return;
1019 }
1020 
1021 
1022 /***********************************************************************//**
1023  * @brief Initialise event bin
1024  *
1025  * This method initialises the event bin. The event bin is cleared and all
1026  * fixed pointers are set.
1027  ***************************************************************************/
1029 {
1030  // Free any existing memory
1031  m_bin.free_members();
1032 
1033  // Set fixed pointers (those will not be set in set_bin)
1034  m_bin.m_time = &m_time;
1035  m_bin.m_ontime = &m_ontime;
1036 
1037  // Return
1038  return;
1039 }
1040 
1041 
1042 /***********************************************************************//**
1043  * @brief Set event bin
1044  *
1045  * @param[in] index Event index [0,...,size()-1].
1046  *
1047  * @exception GException::out_of_range
1048  * Event index is outside valid range.
1049  * @exception GException::invalid_value
1050  * Energy vectors have not been set up.
1051  * Sky directions and solid angles vectors have not been set up.
1052  *
1053  * This method provides the event attributes to the event bin. The event bin
1054  * is in fact physically stored in the event cube, and only a single event
1055  * bin is indeed allocated. This method sets up the pointers in the event
1056  * bin so that a client can easily access the information of individual bins
1057  * as if they were stored in an array.
1058  ***************************************************************************/
1059 void GCTAEventCube::set_bin(const int& index)
1060 {
1061  // Optionally check if the index is valid
1062  #if defined(G_RANGE_CHECK)
1063  if (index < 0 || index >= size()) {
1064  throw GException::out_of_range(G_SET_BIN, "Event index", index, size());
1065  }
1066  #endif
1067 
1068  // Check for the existence of energies and energy widths
1069  if (m_energies.size() != ebins() || m_ewidth.size() != ebins()) {
1070  std::string msg = "Number of energy bins ("+gammalib::str(ebins())+
1071  ") is incompatible with size of energies ("+
1072  gammalib::str(m_energies.size())+") an/or size of "
1073  "energy widths ("+gammalib::str(m_ewidth.size())+
1074  "). Please provide an correctly initialised event "
1075  "cube.";
1077  }
1078 
1079  // Check for the existence of sky directions and solid angles
1080  if (m_dirs.size() != npix() || m_solidangle.size() != npix()) {
1081  std::string msg = "Number of sky pixels ("+gammalib::str(npix())+
1082  ") is incompatible with size of sky directions ("+
1083  gammalib::str(m_dirs.size())+") and/or size of "
1084  "solid angles ("+gammalib::str(m_solidangle.size())+
1085  "). Please provide an correctly initialised event "
1086  "cube.";
1088  }
1089 
1090  // Set pixel and energy bin indices.
1091  m_bin.m_ipix = index % npix();
1092  m_bin.m_ieng = index / npix();
1093 
1094  // Set pointers
1095  m_bin.m_counts = const_cast<double*>(&(m_map.pixels()[index]));
1097  //m_bin.m_time = &m_time;
1098  m_bin.m_dir = &(m_dirs[m_bin.m_ipix]);
1101  //m_bin.m_ontime = &m_ontime;
1102  m_bin.m_weight = const_cast<double*>(&(m_weights.pixels()[index]));
1103 
1104  // Return
1105  return;
1106 }
void clear(void)
Clear CTA pointing.
void read_ebds(const GFitsTable &hdu)
Read energy boundaries from HDU.
const GGti & gti(void) const
Return Good Time Intervals.
Definition: GEvents.hpp:134
Sky map class.
Definition: GSkyMap.hpp:89
Abstract FITS image base class.
Definition: GFitsImage.hpp:43
void init_members(void)
Initialise class members.
Definition: GEventCube.cpp:143
GFitsTable * table(const int &extno)
Get pointer to table HDU.
Definition: GFits.cpp:482
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
const int & ny(void) const
Returns number of pixels in y coordinate.
Definition: GSkyMap.hpp:377
double * m_ontime
Pointer to ontime of bin (seconds)
double * m_counts
Pointer to number of counts.
#define G_SET_DIRECTIONS
virtual GCTAEventBin * operator[](const int &index)
Event bin access operator.
GCTAPointing m_pnt
Event cube pointing.
Definition of GammaLib CTA typemaps.
int size(void) const
Return number of energy boundaries.
Definition: GEbounds.hpp:163
GCTAEventBin class interface definition.
void init_members(void)
Initialise class members.
Definition: GEvents.cpp:176
const GSkyMap & weights(void) const
Return event cube weights as sky map.
virtual GCTAEventCube * clone(void) const
Clone CTA event cube.
void free_members(void)
Delete class members.
Definition: GEventCube.cpp:165
virtual void clear(void)
Clear CTA event cube.
int m_ipix
Index in spatial map.
GCTAEventCube(void)
Void constructor.
virtual GEventCube & operator=(const GEventCube &cube)
Assignment operator.
Definition: GEventCube.cpp:104
void free_members(void)
Delete class members.
const GTime & tstop(void) const
Return stop time.
Definition: GEvents.hpp:158
#define G_SET_ENERGIES
virtual void set_energies(void)
Set log mean energies and energy widths of event cube.
virtual int number(void) const
Return number of events in cube.
void clear(void)
Clear time.
Definition: GTime.cpp:252
GEbounds m_ebounds
Energy boundaries covered by events.
Definition: GEvents.hpp:111
void read(const GFitsTable &table)
Read energy boundaries from FITS table.
Definition: GEbounds.cpp:761
Gammalib tools definition.
FITS file class.
Definition: GFits.hpp:63
const std::string extname_gti
Definition: GGti.hpp:44
FITS file class interface definition.
std::string print(const GChatter &chatter=NORMAL) const
Print sky direction information.
Definition: GSkyDir.cpp:1227
int m_ieng
Index of energy layer.
virtual GCTAEventCube & operator=(const GCTAEventCube &cube)
Assignment operator.
void read(const GFitsHDU &hdu)
Read skymap from FITS HDU.
Definition: GSkyMap.cpp:2664
void write(GFits &file, const std::string &extname=gammalib::extname_ebounds, const std::string &unit="keV") const
Write energy boundaries into FITS object.
Definition: GEbounds.cpp:816
GCTAInstDir * m_dir
Pointer to bin direction.
int size(void) const
Return number of Good Time Intervals.
Definition: GGti.hpp:154
void set_directions(void)
Set sky directions and solid angles of events cube.
int ny(void) const
Return number of bins in Y direction.
GTime m_time
Event cube mean time.
virtual void set_times(void)
Set mean event time and ontime of event cube.
#define G_ENERGY
void saveto(const GFilename &filename, const bool &clobber=false)
Saves to specified FITS file.
Definition: GFits.cpp:1293
std::string print(const GChatter &chatter=NORMAL) const
Print sky map.
Definition: GSkyMap.cpp:2787
std::vector< GEnergy > m_ewidth
Array of energy bin widths.
double real(const std::string &keyname) const
Return card value as double precision.
Definition: GFitsHDU.hpp:423
GTime * m_time
Pointer to bin time.
const std::string extname_cta_weights
Sky map pixel class.
Definition: GSkyPixel.hpp:74
void reference(const GTimeReference &ref)
Set time reference for Good Time Intervals.
Definition: GGti.hpp:257
virtual std::string print(const GChatter &chatter=NORMAL) const
Print event cube information.
CTA event bin container class.
GSkyDir pix2dir(const GSkyPixel &pixel) const
Returns sky direction of pixel.
Definition: GSkyMap.cpp:1360
std::vector< GCTAInstDir > m_dirs
Array of event directions.
int npix(void) const
Return number of pixels in one energy bins of the event cube.
Energy boundaries container class.
Definition: GEbounds.hpp:60
void remove(const int &extno)
Remove HDU from FITS file.
Definition: GFits.cpp:862
double m_ontime
Event cube ontime (sec)
const int & nmaps(void) const
Returns number of maps.
Definition: GSkyMap.hpp:389
GFitsHDU * write(GFits &file, const std::string &extname="") const
Write sky map into FITS file.
Definition: GSkyMap.cpp:2697
Filename class.
Definition: GFilename.hpp:62
void free_members(void)
Delete class members.
Definition: GEvents.cpp:206
int nx(void) const
Return number of bins in X direction.
virtual void load(const GFilename &filename)
Load CTA event cube from FITS file.
GCTAInstDir instdir(const GSkyDir &skydir) const
Get instrument direction from sky direction.
GCTAEventBin m_bin
Actual event bin.
void init_bin(void)
Initialise event bin.
void set_detxy(const GCTAPointing &pnt)
Set DETX and DETY coordinates.
const std::string extname_ebounds
Definition: GEbounds.hpp:44
const std::string extname_cta_counts
CTA pointing class.
Abstract interface for FITS table.
Definition: GFitsTable.hpp:44
GGti m_gti
Good time intervals covered by events.
Definition: GEvents.hpp:112
GChatter
Definition: GTypemaps.hpp:33
void free_members(void)
Delete class members.
const GEnergy & energy(const int &index) const
Return energy of cube layer.
GSkyMap m_weights
Cube weights stored as sky map.
void clear(void)
Clear instance.
Definition: GSkyMap.cpp:1100
Good Time Interval class.
Definition: GGti.hpp:62
void radec_deg(const double &ra, const double &dec)
Set equatorial sky direction (degrees)
Definition: GSkyDir.cpp:224
const GTime & tstop(void) const
Returns latest stop time in Good Time Intervals.
Definition: GGti.hpp:207
const GTime & tstart(void) const
Returns earliest start time in Good Time Intervals.
Definition: GGti.hpp:194
#define G_SET_BIN
std::vector< GEnergy > m_energies
Array of log mean energies.
virtual void clear(void)
Clear eventbin.
void copy_members(const GCTAEventCube &cube)
Copy class members.
double * m_weight
Pointer to weight of bin.
void init_members(void)
Initialise class members.
GFitsImage * image(const int &extno)
Get pointer to image HDU.
Definition: GFits.cpp:368
void read_gti(const GFitsTable &hdu)
Read GTIs from HDU.
const GEnergy & emin(void) const
Return minimum energy.
Definition: GEvents.hpp:170
virtual ~GCTAEventCube(void)
Destructor.
virtual int dim(void) const
Return dimension of event cube.
void read(const GFitsTable &table)
Read Good Time Intervals and time reference from FITS table.
Definition: GGti.cpp:753
GChatter reduce(const GChatter &chatter)
Reduce chattiness by one level.
Definition: GTypemaps.hpp:65
CTA event bin container class interface definition.
Exception handler interface definition.
double * m_solidangle
Pointer to solid angle of pixel (sr)
Implements a time reference.
GEnergy * m_ewidth
Pointer to energy width of bin.
const int & nx(void) const
Returns number of pixels in x coordinate.
Definition: GSkyMap.hpp:361
const GEbounds & ebounds(void) const
Return energy boundaries.
Definition: GEvents.hpp:122
GSkyMap m_map
Counts cube stored as sky map.
CTA instrument direction class.
Definition: GCTAInstDir.hpp:63
const GSkyDir & dir(void) const
Return pointing sky direction.
GEnergy * m_energy
Pointer to bin energy.
#define G_NAXIS
void set_bin(const int &index)
Set event bin.
double solidangle(const int &index) const
Returns solid angle of pixel.
Definition: GSkyMap.cpp:1858
virtual void write(GFits &file) const
Write CTA event cube into FITS file.
#define G_CTA_MJDREF
Reference of CTA time frame.
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
#define G_SET_TIMES
const int & npix(void) const
Returns number of pixels.
Definition: GSkyMap.hpp:345
int ebins(void) const
Return number of energy bins in the event cube.
Sky direction class.
Definition: GSkyDir.hpp:62
void write(GFits &fits, const std::string &extname=gammalib::extname_gti) const
Write Good Time Intervals and time reference into FITS object.
Definition: GGti.cpp:796
Abstract event bin container class.
Definition: GEventCube.hpp:46
void close(void)
Close FITS file.
Definition: GFits.cpp:1342
const GEnergy & emax(void) const
Return maximum energy.
Definition: GEvents.hpp:182
const double & ontime(void) const
Returns ontime.
Definition: GGti.hpp:240
const double * pixels(void) const
Returns pointer to pixel data.
Definition: GSkyMap.hpp:475
virtual void read(const GFits &file)
Read CTA event cube from FITS file.
void read_cntmap(const GFitsImage &hdu)
Read CTA counts map from HDU.
std::vector< double > m_solidangle
Array of solid angles (sr)
GEnergy elogmean(const GEnergy &a, const GEnergy &b)
Computes log mean energy.
Definition: GTools.cpp:1290
Class that handles energies in a unit independent way.
Definition: GEnergy.hpp:48
const GSkyMap & counts(void) const
Return event cube counts as sky map.
bool m_has_pnt
Event cube has pointing.
virtual void save(const GFilename &filename, const bool &clobber=false) const
Save CTA event cube into FITS file.
virtual int naxis(const int &axis) const
Return number of bins in axis.
std::string str(const unsigned short int &value)
Convert unsigned short integer value into string.
Definition: GTools.cpp:489
virtual int size(void) const
Return number of bins in event cube.