GammaLib  1.7.0.dev
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
GLATEventCube.cpp
Go to the documentation of this file.
1 /***************************************************************************
2  * GLATEventCube.cpp - Fermi/LAT event cube class *
3  * ----------------------------------------------------------------------- *
4  * copyright (C) 2009-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 GLATEventCube.cpp
23  * @brief Fermi/LAT event cube 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 "GFilename.hpp"
33 #include "GFitsImage.hpp"
34 #include "GFitsTable.hpp"
35 #include "GLATEventCube.hpp"
36 #include "GLATException.hpp"
37 
38 /* __ Method name definitions ____________________________________________ */
39 #define G_NAXIS "GLATEventCube::naxis(int)"
40 #define G_DIFFNAME "GLATEventCube::diffname(int&)"
41 #define G_DIFFRSP "GLATEventCube::diffrsp(int&)"
42 #define G_READ_SRCMAP "GLATEventCube::read_srcmap(GFitsImage&)"
43 #define G_SET_DIRECTIONS "GLATEventCube::set_directions()"
44 #define G_SET_ENERGIES "GLATEventCube::set_energies()"
45 #define G_SET_TIMES "GLATEventCube::set_times()"
46 #define G_SET_BIN "GLATEventCube::set_bin(int&)"
47 
48 /* __ Macros _____________________________________________________________ */
49 
50 /* __ Coding definitions _________________________________________________ */
51 #define G_MISSING_WCS_KLUDGE //!< Handle source maps with missing WCS
52 
53 /* __ Debug definitions __________________________________________________ */
54 
55 
56 /*==========================================================================
57  = =
58  = Constructors/destructors =
59  = =
60  ==========================================================================*/
61 
62 /***********************************************************************//**
63  * @brief Void constructor
64  ***************************************************************************/
66 {
67  // Initialise members
68  init_members();
69 
70  // Return
71  return;
72 }
73 
74 
75 /***********************************************************************//**
76  * @brief File name constructor
77  *
78  * @param[in] filename Counts cube filename.
79  *
80  * Construct event cube object by loading the events from a FITS file.
81  ***************************************************************************/
83 {
84  // Initialise members
85  init_members();
86 
87  // Load counts cube
88  load(filename);
89 
90  // Return
91  return;
92 }
93 
94 
95 /***********************************************************************//**
96  * @brief Copy constructor
97  *
98  * @param[in] cube LAT event cube.
99  ***************************************************************************/
101 {
102  // Initialise members
103  init_members();
104 
105  // Copy members
106  copy_members(cube);
107 
108  // Return
109  return;
110 }
111 
112 
113 /***********************************************************************//**
114  * @brief Destructor
115  ***************************************************************************/
117 {
118  // Free members
119  free_members();
120 
121  // Return
122  return;
123 }
124 
125 
126 /*==========================================================================
127  = =
128  = Operators =
129  = =
130  ==========================================================================*/
131 
132 /***********************************************************************//**
133  * @brief Assignment operator
134  *
135  * @param[in] cube LAT event cube.
136  * @return LAT event cube.
137  ***************************************************************************/
139 {
140  // Execute only if object is not identical
141  if (this != &cube) {
142 
143  // Copy base class members
144  this->GEventCube::operator=(cube);
145 
146  // Free members
147  free_members();
148 
149  // Initialise members
150  init_members();
151 
152  // Copy members
153  copy_members(cube);
154 
155  } // endif: object was not identical
156 
157  // Return this object
158  return *this;
159 }
160 
161 
162 /***********************************************************************//**
163  * @brief Event bin access operator
164  *
165  * @param[in] index Event index [0,...,size()-1].
166  * @return Pointer to event bin.
167  *
168  * Returns pointer to an event bin.
169  ***************************************************************************/
171 {
172  // Set event bin
173  set_bin(index);
174 
175  // Return pointer
176  return (&m_bin);
177 }
178 
179 
180 /***********************************************************************//**
181  * @brief Event bin access operator (const version)
182  *
183  * @param[in] index Event index [0,...,size()-1].
184  * @return Pointer to event bin.
185  *
186  * Returns pointer to an event bin.
187  ***************************************************************************/
188 const GLATEventBin* GLATEventCube::operator[](const int& index) const
189 {
190  // Set event bin (circumvent const correctness)
191  const_cast<GLATEventCube*>(this)->set_bin(index);
192 
193  // Return pointer
194  return (&m_bin);
195 }
196 
197 
198 /*==========================================================================
199  = =
200  = Public methods =
201  = =
202  ==========================================================================*/
203 
204 /***********************************************************************//**
205  * @brief Clear instance
206  *
207  * This method properly resets the object to an initial state.
208  ***************************************************************************/
210 {
211  // Free class members (base and derived classes, derived class first)
212  free_members();
213  this->GEventCube::free_members();
214  this->GEvents::free_members();
215 
216  // Initialise members
217  this->GEvents::init_members();
218  this->GEventCube::init_members();
219  init_members();
220 
221  // Return
222  return;
223 }
224 
225 
226 /***********************************************************************//**
227  * @brief Clone instance
228  *
229  * @return Deep copy of LAT event cube.
230  ***************************************************************************/
232 {
233  return new GLATEventCube(*this);
234 }
235 
236 
237 /***********************************************************************//**
238  * @brief Return number of bins in event cube
239  *
240  * @return Number of bins in event cube.
241  ***************************************************************************/
242 int GLATEventCube::size(void) const
243 {
244  // Compute number of bins
245  int nbins = m_map.npix() * m_map.nmaps();
246 
247  // Return number of bins
248  return nbins;
249 }
250 
251 
252 /***********************************************************************//**
253  * @brief Return dimension of event cube
254  *
255  * @return Number of dimensions in event cube (either 2 or 3).
256  ***************************************************************************/
257 int GLATEventCube::dim(void) const
258 {
259  // Compute dimension from sky map
260  int dim = (m_map.nmaps() > 1) ? 3 : 2;
261 
262  // Return dimension
263  return dim;
264 }
265 
266 
267 /***********************************************************************//**
268  * @brief Return number of bins in axis
269  *
270  * @param[in] axis Axis [0,...,dim()-1]
271  * @return Number of bins in specified axis
272  *
273  * @exception GException::out_of_range
274  * Axis is out of range.
275  *
276  * Returns the number of bins along a given event cube @p axis.
277  ***************************************************************************/
278 int GLATEventCube::naxis(const int& axis) const
279 {
280  // Optionally check if the axis is valid
281  #if defined(G_RANGE_CHECK)
282  if (axis < 0 || axis >= dim()) {
283  throw GException::out_of_range(G_NAXIS, "LAT event cube axis", axis, dim());
284  }
285  #endif
286 
287  // Set result
288  int naxis = 0;
289  switch (axis) {
290  case 0:
291  naxis = m_map.nx();
292  break;
293  case 1:
294  naxis = m_map.ny();
295  break;
296  case 2:
297  naxis = m_map.npix();
298  break;
299  }
300 
301  // Return result
302  return naxis;
303 }
304 
305 
306 /***********************************************************************//**
307  * @brief Load LAT event cube from FITS file
308  *
309  * @param[in] filename FITS file name.
310  ***************************************************************************/
311 void GLATEventCube::load(const GFilename& filename)
312 {
313  // Open FITS file
314  GFits fits(filename);
315 
316  // Read event cube from FITS file
317  read(fits);
318 
319  // Close FITS file
320  fits.close();
321 
322  // Return
323  return;
324 }
325 
326 
327 /***********************************************************************//**
328  * @brief Save LAT event cube into FITS file
329  *
330  * @param[in] filename FITS file name.
331  * @param[in] clobber Overwrite existing FITS file? (default: false)
332  *
333  * Save the LAT event cube into FITS file.
334  ***************************************************************************/
335 void GLATEventCube::save(const GFilename& filename,
336  const bool& clobber) const
337 {
338  // Create empty FITS file
339  GFits fits;
340 
341  // Write event cube into FITS file
342  write(fits);
343 
344  // Save FITS file
345  fits.saveto(filename, clobber);
346 
347  // Return
348  return;
349 }
350 
351 
352 /***********************************************************************//**
353  * @brief Read LAT event cube from FITS file.
354  *
355  * @param[in] fits FITS file.
356  *
357  * It is assumed that the counts map resides in the primary extension of the
358  * FITS file, the energy boundaries reside in the `EBOUNDS` extension and the
359  * Good Time Intervals reside in the `GTI` extension. The method clears the
360  * object before loading, thus any events residing in the object before
361  * loading will be lost.
362  ***************************************************************************/
363 void GLATEventCube::read(const GFits& fits)
364 {
365  // Clear object
366  clear();
367 
368  // Get HDUs
369  const GFitsImage& hdu_cntmap = *fits.image("Primary");
370  const GFitsTable& hdu_ebounds = *fits.table(gammalib::extname_ebounds);
371  const GFitsTable& hdu_gti = *fits.table(gammalib::extname_gti);
372 
373  // Load counts map
374  read_cntmap(hdu_cntmap);
375 
376  // Load energy boundaries
377  read_ebds(hdu_ebounds);
378 
379  // Load GTIs
380  read_gti(hdu_gti);
381 
382  // Load additional source maps
383  for (int i = 1; i < fits.size(); ++i) {
384 
385  // Only consider image extensions
386  if (fits.at(i)->exttype() == GFitsHDU::HT_IMAGE) {
387 
388  // Get FITS image
389  const GFitsImage& hdu_srcmap = *fits.image(i);
390 
391  // Handle files that do not have WCS in their additional
392  // source map extensions
393  #if defined(G_MISSING_WCS_KLUDGE)
394  std::string keys[] = {"CTYPE1", "CTYPE2", "CTYPE3",
395  "CRPIX1", "CRPIX2", "CRPIX3",
396  "CRVAL1", "CRVAL2", "CRVAL3",
397  "CDELT1", "CDELT2", "CDELT3",
398  "CUNIT1", "CUNIT2", "CUNIT3",
399  "CROTA2"};
400  for (int i = 0; i < 16; ++i) {
401  if (!hdu_srcmap.has_card(keys[i])) {
402  const_cast<GFitsImage&>(hdu_srcmap).card(hdu_cntmap.card(keys[i]));
403  }
404  }
405  #endif
406 
407  // Read source map
408  read_srcmap(hdu_srcmap);
409 
410  } // endif: extension is an image
411 
412  } // endfor: looped over additional source maps
413 
414  // Return
415  return;
416 }
417 
418 
419 /***********************************************************************//**
420  * @brief Write LAT event cube into FITS file
421  *
422  * @param[in] fits FITS file.
423  ***************************************************************************/
424 void GLATEventCube::write(GFits& fits) const
425 {
426  // Write cube
427  m_map.write(fits);
428 
429  // Write energy boundaries
430  ebounds().write(fits);
431 
432  // Write Good Time intervals
433  gti().write(fits);
434 
435  // Write additional source maps
436  for (int i = 0; i < m_srcmap.size(); ++i) {
437  m_srcmap[i]->write(fits);
438  }
439 
440  // Return
441  return;
442 }
443 
444 
445 /***********************************************************************//**
446  * @brief Return number of events in cube
447  *
448  * @return Total number of events in event cube.
449  ***************************************************************************/
450 int GLATEventCube::number(void) const
451 {
452  // Initialise result
453  double number = 0.0;
454 
455  // Get pointer on skymap pixels
456  const double* pixels = m_map.pixels();
457 
458  // Sum event cube
459  if (size() > 0 && pixels != NULL) {
460  for (int i = 0; i < size(); ++i) {
461  number += pixels[i];
462  }
463  }
464 
465  // Return
466  return int(number+0.5);
467 }
468 
469 
470 /***********************************************************************//**
471  * @brief Set event cube from sky map
472  ***************************************************************************/
473 void GLATEventCube::map(const GSkyMap& map)
474 {
475  // Store sky map
476  m_map = map;
477 
478  // Compute sky directions
479  set_directions();
480 
481  // Return
482  return;
483 }
484 
485 
486 /***********************************************************************//**
487  * @brief Print event cube information
488  *
489  * @param[in] chatter Chattiness (defaults to NORMAL).
490  * @return String containing event cube information.
491  ***************************************************************************/
492 std::string GLATEventCube::print(const GChatter& chatter) const
493 {
494  // Initialise result string
495  std::string result;
496 
497  // Continue only if chatter is not silent
498  if (chatter != SILENT) {
499 
500  // Append header
501  result.append("=== GLATEventCube ===");
502 
503  // Append information
504  result.append("\n"+gammalib::parformat("Number of elements") +
505  gammalib::str(size()));
506  result.append("\n"+gammalib::parformat("Number of pixels"));
507  result.append(gammalib::str(m_map.nx()) +
508  " x " +
509  gammalib::str(m_map.ny()));
510  result.append("\n"+gammalib::parformat("Number of energy bins") +
511  gammalib::str(ebins()));
512  result.append("\n"+gammalib::parformat("Number of events") +
513  gammalib::str(number()));
514 
515  // Append time interval
516  result.append("\n"+gammalib::parformat("Time interval"));
517  if (gti().size() > 0) {
518  result.append(gammalib::str(tstart().secs()) +
519  " - " +
520  gammalib::str(tstop().secs())+" sec");
521  }
522  else {
523  result.append("not defined");
524  }
525 
526  // Append energy range
527  result.append("\n"+gammalib::parformat("Energy range"));
528  if (ebounds().size() > 0) {
529  result.append(emin().print()+" - "+emax().print(chatter));
530  }
531  else {
532  result.append("not defined");
533  }
534 
535  // Append detailed information
536  GChatter reduced_chatter = gammalib::reduce(chatter);
537  if (reduced_chatter > SILENT) {
538 
539  // Append sky projection
540  if (m_map.projection() != NULL) {
541  result.append("\n"+m_map.projection()->print(reduced_chatter));
542  }
543 
544  // Append source maps
545  result.append("\n"+gammalib::parformat("Number of source maps"));
546  result.append(gammalib::str(m_srcmap.size()));
547  for (int i = 0; i < m_srcmap.size(); ++i) {
548  result.append("\n"+gammalib::parformat(" "+m_srcmap_names[i]));
549  result.append(gammalib::str(m_srcmap[i]->nx()));
550  result.append(" x ");
551  result.append(gammalib::str(m_srcmap[i]->ny()));
552  result.append(" x ");
553  result.append(gammalib::str(m_srcmap[i]->nmaps()));
554  }
555 
556  }
557 
558  } // endif: chatter was not silent
559 
560  // Return result
561  return result;
562 }
563 
564 
565 /***********************************************************************//**
566  * @brief Return name of diffuse model
567  *
568  * @param[in] index Diffuse model index [0,...,ndiffrsp()-1].
569  * @return Name of diffuse model.
570  *
571  * @exception GException::out_of_range
572  * Model index out of valid range.
573  *
574  * Returns name of diffuse model.
575  ***************************************************************************/
576 std::string GLATEventCube::diffname(const int& index) const
577 {
578  // Optionally check if the index is valid
579  #if defined(G_RANGE_CHECK)
580  if (index < 0 || index >= ndiffrsp()) {
581  throw GException::out_of_range(G_DIFFNAME, index, 0, ndiffrsp()-1);
582  }
583  #endif
584 
585  // Return
586  return m_srcmap_names[index];
587 }
588 
589 
590 /***********************************************************************//**
591  * @brief Return diffuse response map
592  *
593  * @param[in] index Diffuse model index [0,...,ndiffrsp()-1].
594  * @return Pointer to diffuse response map.
595  *
596  * @exception GException::out_of_range
597  * Model index out of valid range.
598  *
599  * Returns pointer to diffuse model sky map.
600  ***************************************************************************/
601 GSkyMap* GLATEventCube::diffrsp(const int& index) const
602 {
603  // Optionally check if the index is valid
604  #if defined(G_RANGE_CHECK)
605  if (index < 0 || index >= ndiffrsp()) {
606  throw GException::out_of_range(G_DIFFRSP, index, 0, ndiffrsp()-1);
607  }
608  #endif
609 
610  // Return
611  return m_srcmap[index];
612 }
613 
614 
615 /***********************************************************************//**
616  * @brief Computes the maximum radius (in degrees) around a given source
617  * direction that fits spatially into the event cube
618  *
619  * @param[in] srcDir Source direction.
620  * @return Maximum radius in degrees that fully fits into event cube.
621  *
622  * By computing the sky directions of the event cube boundaries, the maximum
623  * radius is computed that fits fully within the event cube. This method is
624  * used for PSF normalization.
625  ***************************************************************************/
626 double GLATEventCube::maxrad(const GSkyDir& srcDir) const
627 {
628  // Initialise radius
629  double radius = 0.0;
630 
631  // Continue only if sky direction is within sky map
632  if (m_map.contains(srcDir)) {
633 
634  // Set to largest possible radius
635  radius = 180.0;
636 
637  // Move along upper edge in longitude
638  int iy = 0;
639  for (int ix = 0; ix < nx(); ++ix) {
640  GSkyPixel pixel = GSkyPixel(double(ix), double(iy));
641  double distance = m_map.pix2dir(pixel).dist_deg(srcDir);
642  if (distance < radius) {
643  radius = distance;
644  }
645  }
646 
647  // Move along lower edge in longitude
648  iy = ny()-1;
649  for (int ix = 0; ix < nx(); ++ix) {
650  GSkyPixel pixel = GSkyPixel(double(ix), double(iy));
651  double distance = m_map.pix2dir(pixel).dist_deg(srcDir);
652  if (distance < radius) {
653  radius = distance;
654  }
655  }
656 
657  // Move along left edge in latitude
658  int ix = 0;
659  for (int iy = 0; iy < ny(); ++iy) {
660  GSkyPixel pixel = GSkyPixel(double(ix), double(iy));
661  double distance = m_map.pix2dir(pixel).dist_deg(srcDir);
662  if (distance < radius) {
663  radius = distance;
664  }
665  }
666 
667  // Move along right edge in latitude
668  ix = nx()-1;
669  for (int iy = 0; iy < ny(); ++iy) {
670  GSkyPixel pixel = GSkyPixel(double(ix), double(iy));
671  double distance = m_map.pix2dir(pixel).dist_deg(srcDir);
672  if (distance < radius) {
673  radius = distance;
674  }
675  }
676 
677  } // endif: sky direction within sky map
678 
679  // Return radius
680  return radius;
681 }
682 
683 
684 /*==========================================================================
685  = =
686  = Private methods =
687  = =
688  ==========================================================================*/
689 
690 /***********************************************************************//**
691  * @brief Initialise class members
692  ***************************************************************************/
694 {
695  // Initialise members
696  m_bin.clear();
697  m_map.clear();
698  m_time.clear();
699  m_srcmap.clear();
700  m_srcmap_names.clear();
701  m_enodes.clear();
702  m_dirs.clear();
703  m_solidangle.clear();
704  m_energies.clear();
705  m_ewidth.clear();
706  m_ontime = 0.0;
707 
708  // Return
709  return;
710 }
711 
712 
713 /***********************************************************************//**
714  * @brief Copy class members
715  *
716  * @param[in] cube Fermi/LAT event cube.
717  ***************************************************************************/
719 {
720  // Copy Fermi/LAT specific attributes
721  m_bin = cube.m_bin;
722  m_map = cube.m_map;
723  m_time = cube.m_time;
724  m_ontime = cube.m_ontime;
725  m_enodes = cube.m_enodes;
726  m_dirs = cube.m_dirs;
727  m_solidangle = cube.m_solidangle;
728  m_energies = cube.m_energies;
729  m_ewidth = cube.m_ewidth;
730 
731  // Clone source maps and copy source map names
732  m_srcmap.clear();
733  for (int i = 0; i < cube.m_srcmap.size(); ++i) {
734  m_srcmap.push_back(cube.m_srcmap[i]->clone());
735  }
737 
738  // Return
739  return;
740 }
741 
742 
743 /***********************************************************************//**
744  * @brief Delete class members
745  ***************************************************************************/
747 {
748  // Release source maps
749  for (int i = 0; i < m_srcmap.size(); ++i) {
750  if (m_srcmap[i] != NULL) delete m_srcmap[i];
751  m_srcmap[i] = NULL;
752  }
753  m_srcmap.clear();
754  m_srcmap_names.clear();
755 
756  // Return
757  return;
758 }
759 
760 
761 /***********************************************************************//**
762  * @brief Read Fermi/LAT counts map from HDU.
763  *
764  * @param[in] hdu Image HDU.
765  *
766  * This method reads a Fermi/LAT counts map from a FITS image. The counts map
767  * is stored in a GSkyMap object, and a pointer is set up to access the
768  * pixels individually. Recall that skymap pixels are stored in the order
769  * (ix,iy,ebin).
770  ***************************************************************************/
772 {
773  // Load counts map as sky map
774  m_map.read(hdu);
775 
776  // Set sky directions
777  set_directions();
778 
779  // Return
780  return;
781 }
782 
783 
784 /***********************************************************************//**
785  * @brief Read LAT source map from HDU.
786  *
787  * @param[in] hdu Image HDU.
788  *
789  * @exception GLATException::wcs_incompatible
790  * Source map not compatible with sky map
791  *
792  * This method reads a LAT source map from a FITS image. The source map is
793  * stored in a GSkyMap object and is given in units of counts/pixel/MeV.
794  ***************************************************************************/
796 {
797  // Allocate skymap
798  GSkyMap* map = new GSkyMap;
799 
800  // Read skymap
801  map->read(hdu);
802 
803  // Check that source map sky projection is consistent with counts
804  // map sky projection
805  if (*(m_map.projection()) != *(map->projection())) {
807  }
808 
809  // Check that source map dimension is consistent with counts map
810  // dimension
811  if (m_map.nx() != map->nx() ||
812  m_map.ny() != map->ny()) {
814  }
815 
816  // Check that source map has required number of energy bins
817  if (m_map.nmaps()+1 != map->nmaps()) {
819  }
820 
821  // Append source map to list of maps
822  m_srcmap.push_back(map);
823  m_srcmap_names.push_back(hdu.extname());
824 
825  // Return
826  return;
827 }
828 
829 
830 /***********************************************************************//**
831  * @brief Read energy boundaries from HDU.
832  *
833  * @param[in] hdu Energy boundaries table.
834  *
835  * Read the energy boundaries from the HDU.
836  *
837  * @todo Energy bounds read method should take const GFitsTable* as argument
838  ***************************************************************************/
840 {
841  // Read energy boundaries
842  m_ebounds.read(hdu);
843 
844  // Set log mean energies and energy widths
845  set_energies();
846 
847  // Return
848  return;
849 }
850 
851 
852 /***********************************************************************//**
853  * @brief Read GTIs from HDU.
854  *
855  * @param[in] hdu GTI table.
856  *
857  * Reads the Good Time Intervals from the GTI extension. Since the Fermi
858  * LAT Science Tools do not set corrently the time reference for source
859  * maps, the method automatically adds this missing information so that
860  * the time reference is set correctly. The time reference that is assumed
861  * for Fermi LAT is
862  *
863  * MJDREFI 51910
864  * MJDREFF 0.00074287037037037
865  *
866  ***************************************************************************/
868 {
869  // Work on a local copy of the HDU to make the kluge work
870  GFitsTable* hdu_local = hdu.clone();
871 
872  // Kluge: modify HDU table in case that the MJDREF header keyword is
873  // blank. This happens for Fermi LAT source maps since the Science
874  // Tools do not properly write out the time reference. We hard-code
875  // here the Fermi LAT time reference to circumvent the problem.
876  if (hdu.has_card("MJDREF")) {
877  if (gammalib::strip_whitespace(hdu.string("MJDREF")).empty()) {
878  const_cast<GFitsHeader&>(hdu_local->header()).remove("MJDREF");
879  hdu_local->card("MJDREFI", 51910,
880  "Integer part of MJD reference");
881  hdu_local->card("MJDREFF", 0.00074287037037037,
882  "Fractional part of MJD reference");
883  }
884  }
885 
886  // Read Good Time Intervals
887  m_gti.read(*hdu_local);
888 
889  // Set time
890  set_times();
891 
892  // Return
893  return;
894 }
895 
896 
897 /***********************************************************************//**
898  * @brief Set sky directions and solid angles of events cube.
899  *
900  * @exception GLATException::no_sky
901  * No sky pixels found in event cube.
902  *
903  * This method computes the sky directions and solid angles for all event
904  * cube pixels. Sky directions are stored in an array of GLATInstDir objects
905  * while solid angles are stored in units of sr in a double precision array.
906  ***************************************************************************/
908 {
909  // Throw an error if we have no sky pixels
910  if (npix() < 1) {
911  throw GLATException::no_sky(G_SET_DIRECTIONS, "Every LAT event cube"
912  " needs a definition of the sky pixels.");
913  }
914 
915  // Clear old pixel directions and solid angle
916  m_dirs.clear();
917  m_solidangle.clear();
918 
919  // Reserve space for pixel directions and solid angles
920  m_dirs.reserve(npix());
921  m_solidangle.reserve(npix());
922 
923  // Set pixel directions and solid angles
924  for (int iy = 0; iy < ny(); ++iy) {
925  for (int ix = 0; ix < nx(); ++ix) {
926  GSkyPixel pixel = GSkyPixel(double(ix), double(iy));
927  m_dirs.push_back(GLATInstDir(m_map.pix2dir(pixel)));
928  m_solidangle.push_back(m_map.solidangle(pixel));
929  }
930  }
931 
932  // Return
933  return;
934 }
935 
936 
937 /***********************************************************************//**
938  * @brief Set log mean energies and energy widths of event cube.
939  *
940  * @exception GLATException::no_ebds
941  * No energy boundaries found in event cube.
942  *
943  * This method computes the log mean energies and the energy widths of the
944  * event cube. The log mean energies and energy widths are stored unit
945  * independent in arrays of GEnergy objects.
946  ***************************************************************************/
948 {
949  // Throw an error if we have no energy bins
950  if (ebins() < 1) {
951  throw GLATException::no_ebds(G_SET_ENERGIES, "Every LAT event cube"
952  " needs a definition of the energy boundaries.");
953  }
954 
955  // Clear old bin energies and energy widths
956  m_energies.clear();
957  m_ewidth.clear();
958  m_enodes.clear();
959 
960  // Reserve space for bin energies and energy widths
961  m_energies.reserve(ebins());
962  m_ewidth.reserve(ebins());
963 
964  // Setup bin energies, energy widths and energy nodes
965  for (int i = 0; i < ebins(); ++i) {
966  m_energies.push_back(ebounds().elogmean(i));
967  m_ewidth.push_back(ebounds().emax(i) - ebounds().emin(i));
968  m_enodes.append(log10(ebounds().emin(i).MeV()));
969  }
970  m_enodes.append(log10(ebounds().emax(ebins()-1).MeV()));
971 
972  // Return
973  return;
974 }
975 
976 
977 /***********************************************************************//**
978  * @brief Set mean event time and ontime of event cube
979  *
980  * @exception GException::invalid_value
981  * No Good Time Intervals found.
982  *
983  * Computes the mean time of the event cube by taking the mean between start
984  * and stop time. Computes also the ontime by summing up of all good time
985  * intervals.
986  *
987  * @todo Could add a more sophisticated mean event time computation that
988  * weights by the length of the GTIs, yet so far we do not really use
989  * the mean event time, hence there is no rush to implement this.
990  ***************************************************************************/
992 {
993  // Throw an error if GTI is empty
994  if (m_gti.size() < 1) {
995  std::string msg = "No Good Time Intervals have been found in event "
996  "cube. Every LAT event cube needs a definition "
997  "of the Good Time Intervals.";
999  }
1000 
1001  // Compute mean time
1002  m_time = m_gti.tstart() + 0.5 * (m_gti.tstop() - m_gti.tstart());
1003 
1004  // Set ontime
1005  m_ontime = m_gti.ontime();
1006 
1007  // Return
1008  return;
1009 }
1010 
1011 
1012 /***********************************************************************//**
1013  * @brief Set event bin
1014  *
1015  * @param[in] index Event index [0,...,size()-1].
1016  *
1017  * @exception GException::out_of_range
1018  * Event index is outside valid range.
1019  * @exception GLATException::no_energies
1020  * Energy vectors have not been set up.
1021  * @exception GLATException::no_dirs
1022  * Sky directions and solid angles vectors have not been set up.
1023  *
1024  * This method provides the event attributes to the event bin. The event bin
1025  * is in fact physically stored in the event cube, and only a single event
1026  * bin is indeed allocated. This method sets up the pointers in the event
1027  * bin so that a client can easily access the information of individual bins
1028  * as if they were stored in an array.
1029  ***************************************************************************/
1030 void GLATEventCube::set_bin(const int& index)
1031 {
1032  // Optionally check if the index is valid
1033  #if defined(G_RANGE_CHECK)
1034  if (index < 0 || index >= size()) {
1035  throw GException::out_of_range(G_SET_BIN, index, 0, size()-1);
1036  }
1037  #endif
1038 
1039  // Check for the existence of energies and energy widths
1040  if (m_energies.size() != ebins() || m_ewidth.size() != ebins()) {
1042  }
1043 
1044  // Check for the existence of sky directions and solid angles
1045  if (m_dirs.size() != npix() || m_solidangle.size() != npix()) {
1047  }
1048 
1049  // Get pixel and energy bin indices.
1050  m_bin.m_index = index;
1051  m_bin.m_ipix = index % npix();
1052  m_bin.m_ieng = index / npix();
1053 
1054  // Set pointers
1055  m_bin.m_cube = this;
1056  m_bin.m_counts = const_cast<double*>(&(m_map.pixels()[index]));
1058  m_bin.m_time = &m_time;
1059  m_bin.m_dir = &(m_dirs[m_bin.m_ipix]);
1062  m_bin.m_ontime = &m_ontime;
1063 
1064  // Return
1065  return;
1066 }
GTime m_time
Event cube mean time.
GFitsHDU * at(const int &extno)
Get pointer to HDU.
Definition: GFits.cpp:233
virtual int naxis(const int &axis) const
Return number of bins in axis.
const GGti & gti(void) const
Return Good Time Intervals.
Definition: GEvents.hpp:134
Sky map class.
Definition: GSkyMap.hpp:89
double dist_deg(const GSkyDir &dir) const
Compute angular distance between sky directions in degrees.
Definition: GSkyDir.hpp:280
GSkyMap * diffrsp(const int &index) const
Return diffuse response map.
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 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
const GSkyProjection * projection(void) const
Returns pointer to sky projection.
Definition: GSkyMap.hpp:414
virtual void set_times(void)
Set mean event time and ontime of event cube.
virtual ~GLATEventCube(void)
Destructor.
int m_index
Actual skymap index.
int ndiffrsp(void) const
Return number of diffuse model components.
void free_members(void)
Delete class members.
void init_members(void)
Initialise class members.
Definition: GEvents.cpp:176
Fermi/LAT instrument direction class.
Definition: GLATInstDir.hpp:44
double maxrad(const GSkyDir &dir) const
Computes the maximum radius (in degrees) around a given source direction that fits spatially into the...
void free_members(void)
Delete class members.
Definition: GEventCube.cpp:165
virtual void load(const GFilename &filename)
Load LAT event cube from FITS file.
void read_gti(const GFitsTable &hdu)
Read GTIs from HDU.
virtual GEventCube & operator=(const GEventCube &cube)
Assignment operator.
Definition: GEventCube.cpp:104
const GTime & tstop(void) const
Return stop time.
Definition: GEvents.hpp:158
int nx(void) const
Return number of bins in X direction.
void clear(void)
Clear node array.
Definition: GNodeArray.cpp:235
virtual HDUType exttype(void) const =0
void clear(void)
Clear time.
Definition: GTime.cpp:251
GEbounds m_ebounds
Energy boundaries covered by events.
Definition: GEvents.hpp:111
std::vector< GEnergy > m_ewidth
Array of energy bin widths.
void read(const GFitsTable &table)
Read energy boundaries from FITS table.
Definition: GEbounds.cpp:776
Gammalib tools definition.
virtual void read(const GFits &file)
Read LAT event cube from FITS file.
FITS file class.
Definition: GFits.hpp:63
const std::string extname_gti
Definition: GGti.hpp:44
const GFitsHeader & header(void) const
Return extension header.
Definition: GFitsHDU.hpp:205
GLATInstDir * m_dir
Pointer to bin direction.
GEnergy * m_ewidth
Pointer to energy width of bin.
std::string strip_whitespace(const std::string &arg)
Strip leading and trailing whitespace from string.
Definition: GTools.cpp:73
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
int size(void) const
Return number of Good Time Intervals.
Definition: GGti.hpp:153
virtual int size(void) const
Return number of bins in event cube.
Abstract FITS image base class definition.
virtual void clear(void)
Clear instance.
GEnergy * m_energy
Pointer to bin energy.
std::vector< double > m_solidangle
Array of solid angles (sr)
int m_ipix
Actual spatial index.
void read_srcmap(const GFitsImage &hdu)
Read LAT source map from HDU.
void set_bin(const int &index)
Set event bin.
Fermi/LAT event bin class.
void saveto(const GFilename &filename, const bool &clobber=false)
Saves to specified FITS file.
Definition: GFits.cpp:1267
int m_ieng
Actual energy index.
Sky map pixel class.
Definition: GSkyPixel.hpp:74
GNodeArray m_enodes
Energy nodes.
void init_members(void)
Initialise class members.
Interface for FITS header class.
Definition: GFitsHeader.hpp:49
#define G_SET_ENERGIES
GSkyDir pix2dir(const GSkyPixel &pixel) const
Returns sky direction of pixel.
Definition: GSkyMap.cpp:1350
GLATEventCube(void)
Void constructor.
virtual void save(const GFilename &filename, const bool &clobber=false) const
Save LAT event cube into FITS file.
const int & nmaps(void) const
Returns number of maps.
Definition: GSkyMap.hpp:378
GLATEventBin m_bin
Actual energy bin.
int npix(void) const
Return number of pixels in event cube sky map.
double * m_solidangle
Pointer to solid angle of pixel (sr)
int ebins(void) const
Return number of energy bins in event cube.
LAT exception handler interface definition.
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
#define G_NAXIS
virtual GLATEventCube & operator=(const GLATEventCube &cube)
Assignment operator.
const std::string extname_ebounds
Definition: GEbounds.hpp:44
bool contains(const GSkyDir &dir) const
Checks if sky direction falls in map.
Definition: GSkyMap.cpp:1879
virtual void set_energies(void)
Set log mean energies and energy widths of event cube.
Abstract interface for FITS table.
Definition: GFitsTable.hpp:44
GGti m_gti
Good time intervals covered by events.
Definition: GEvents.hpp:112
virtual void write(GFits &file) const
Write LAT event cube into FITS file.
GChatter
Definition: GTypemaps.hpp:33
Fermi/LAT event cube class definition.
void copy_members(const GLATEventCube &cube)
Copy class members.
double m_ontime
Event cube ontime (sec)
void clear(void)
Clear instance.
Definition: GSkyMap.cpp:1090
Fermi/LAT event cube class.
const GTime & tstop(void) const
Returns latest stop time in Good Time Intervals.
Definition: GGti.hpp:206
void read_ebds(const GFitsTable &hdu)
Read energy boundaries from HDU.
const std::string & extname(void) const
Return extension name.
Definition: GFitsHDU.hpp:162
#define G_SET_DIRECTIONS
#define G_DIFFRSP
const GTime & tstart(void) const
Returns earliest start time in Good Time Intervals.
Definition: GGti.hpp:193
virtual std::string print(const GChatter &chatter=NORMAL) const
Print event cube information.
#define G_SET_TIMES
std::vector< GSkyMap * > m_srcmap
Pointers to source maps.
std::vector< GEnergy > m_energies
Array of log mean energies.
virtual std::string print(const GChatter &chatter=NORMAL) const =0
Print content of object.
#define G_SET_BIN
virtual void clear(void)
Clear event bin.
double * m_ontime
Pointer to ontime of bin (seconds)
double * m_counts
Pointer to number of counts.
#define G_READ_SRCMAP
GFitsImage * image(const int &extno)
Get pointer to image HDU.
Definition: GFits.cpp:360
void set_directions(void)
Set sky directions and solid angles of events cube.
const GEnergy & emin(void) const
Return minimum energy.
Definition: GEvents.hpp:170
void read(const GFitsTable &table)
Read Good Time Intervals and time reference from FITS table.
Definition: GGti.cpp:639
int size(void) const
Return number of HDUs in FITS file.
Definition: GFits.hpp:237
GChatter reduce(const GChatter &chatter)
Reduce chattiness by one level.
Definition: GTypemaps.hpp:65
std::vector< GLATInstDir > m_dirs
Array of event directions.
std::string string(const std::string &keyname) const
Return card value as string.
Definition: GFitsHDU.hpp:410
virtual GFitsTable * clone(void) const =0
Clones object.
#define G_DIFFNAME
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
virtual GLATEventBin * operator[](const int &index)
Event bin access operator.
const GSkyMap & map(void) const
Return event cube sky map.
virtual int number(void) const
Return number of events in cube.
virtual GLATEventCube * clone(void) const
Clone instance.
virtual int dim(void) const
Return dimension of event cube.
GLATEventCube * m_cube
Event cube back pointer.
void read_cntmap(const GFitsImage &hdu)
Read Fermi/LAT counts map from HDU.
double solidangle(const int &index) const
Returns solid angle of pixel.
Definition: GSkyMap.cpp:1770
GTime * m_time
Pointer to bin time.
std::string parformat(const std::string &s, const int &indent=0)
Convert string in parameter format.
Definition: GTools.cpp:1022
void append(const double &node)
Append one node to array.
Definition: GNodeArray.cpp:351
const GTime & tstart(void) const
Return start time.
Definition: GEvents.hpp:146
const int & npix(void) const
Returns number of pixels.
Definition: GSkyMap.hpp:334
GFitsHeaderCard & card(const int &cardno)
Return header card.
Definition: GFitsHDU.hpp:259
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
int ny(void) const
Return number of bins in Y direction.
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
std::vector< std::string > m_srcmap_names
Source map names.
Filename class interface definition.
GSkyMap m_map
Counts map stored as sky map.
GVector log10(const GVector &vector)
Computes base10 logarithm of vector elements.
Definition: GVector.cpp:1205
GEnergy elogmean(const GEnergy &a, const GEnergy &b)
Computes log mean energy.
Definition: GTools.cpp:1169
std::string diffname(const int &index) const
Return name of diffuse model.
std::string str(const unsigned short int &value)
Convert unsigned short integer value into string.
Definition: GTools.cpp:413
FITS table abstract base class interface definition.