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