GammaLib  2.0.0
 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-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 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 "GException.hpp"
32 #include "GTools.hpp"
33 #include "GFilename.hpp"
34 #include "GFitsImage.hpp"
35 #include "GFitsTable.hpp"
36 #include "GLATEventCube.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",
284  axis, dim());
285  }
286  #endif
287 
288  // Set result
289  int naxis = 0;
290  switch (axis) {
291  case 0:
292  naxis = m_map.nx();
293  break;
294  case 1:
295  naxis = m_map.ny();
296  break;
297  case 2:
298  naxis = m_map.npix();
299  break;
300  }
301 
302  // Return result
303  return naxis;
304 }
305 
306 
307 /***********************************************************************//**
308  * @brief Load LAT event cube from FITS file
309  *
310  * @param[in] filename FITS file name.
311  ***************************************************************************/
312 void GLATEventCube::load(const GFilename& filename)
313 {
314  // Open FITS file
315  GFits fits(filename);
316 
317  // Read event cube from FITS file
318  read(fits);
319 
320  // Close FITS file
321  fits.close();
322 
323  // Return
324  return;
325 }
326 
327 
328 /***********************************************************************//**
329  * @brief Save LAT event cube into FITS file
330  *
331  * @param[in] filename FITS file name.
332  * @param[in] clobber Overwrite existing FITS file? (default: false)
333  *
334  * Save the LAT event cube into FITS file.
335  ***************************************************************************/
336 void GLATEventCube::save(const GFilename& filename,
337  const bool& clobber) const
338 {
339  // Create empty FITS file
340  GFits fits;
341 
342  // Write event cube into FITS file
343  write(fits);
344 
345  // Save FITS file
346  fits.saveto(filename, clobber);
347 
348  // Return
349  return;
350 }
351 
352 
353 /***********************************************************************//**
354  * @brief Read LAT event cube from FITS file.
355  *
356  * @param[in] fits FITS file.
357  *
358  * It is assumed that the counts map resides in the primary extension of the
359  * FITS file, the energy boundaries reside in the `EBOUNDS` extension and the
360  * Good Time Intervals reside in the `GTI` extension. The method clears the
361  * object before loading, thus any events residing in the object before
362  * loading will be lost.
363  ***************************************************************************/
364 void GLATEventCube::read(const GFits& fits)
365 {
366  // Clear object
367  clear();
368 
369  // Get HDUs
370  const GFitsImage& hdu_cntmap = *fits.image("Primary");
371  const GFitsTable& hdu_ebounds = *fits.table(gammalib::extname_ebounds);
372  const GFitsTable& hdu_gti = *fits.table(gammalib::extname_gti);
373 
374  // Load counts map
375  read_cntmap(hdu_cntmap);
376 
377  // Load energy boundaries
378  read_ebds(hdu_ebounds);
379 
380  // Load GTIs
381  read_gti(hdu_gti);
382 
383  // Load additional source maps
384  for (int i = 1; i < fits.size(); ++i) {
385 
386  // Only consider image extensions
387  if (fits.at(i)->exttype() == GFitsHDU::HT_IMAGE) {
388 
389  // Get FITS image
390  const GFitsImage& hdu_srcmap = *fits.image(i);
391 
392  // Handle files that do not have WCS in their additional
393  // source map extensions
394  #if defined(G_MISSING_WCS_KLUDGE)
395  std::string keys[] = {"CTYPE1", "CTYPE2", "CTYPE3",
396  "CRPIX1", "CRPIX2", "CRPIX3",
397  "CRVAL1", "CRVAL2", "CRVAL3",
398  "CDELT1", "CDELT2", "CDELT3",
399  "CUNIT1", "CUNIT2", "CUNIT3",
400  "CROTA2"};
401  for (int i = 0; i < 16; ++i) {
402  if (!hdu_srcmap.has_card(keys[i])) {
403  const_cast<GFitsImage&>(hdu_srcmap).card(hdu_cntmap.card(keys[i]));
404  }
405  }
406  #endif
407 
408  // Read source map
409  read_srcmap(hdu_srcmap);
410 
411  } // endif: extension is an image
412 
413  } // endfor: looped over additional source maps
414 
415  // Return
416  return;
417 }
418 
419 
420 /***********************************************************************//**
421  * @brief Write LAT event cube into FITS file
422  *
423  * @param[in] fits FITS file.
424  ***************************************************************************/
425 void GLATEventCube::write(GFits& fits) const
426 {
427  // Write cube
428  m_map.write(fits);
429 
430  // Write energy boundaries
431  ebounds().write(fits);
432 
433  // Write Good Time intervals
434  gti().write(fits);
435 
436  // Write additional source maps
437  for (int i = 0; i < m_srcmap.size(); ++i) {
438  m_srcmap[i]->write(fits);
439  }
440 
441  // Return
442  return;
443 }
444 
445 
446 /***********************************************************************//**
447  * @brief Return number of events in cube
448  *
449  * @return Total number of events in event cube.
450  ***************************************************************************/
451 int GLATEventCube::number(void) const
452 {
453  // Initialise result
454  double number = 0.0;
455 
456  // Get pointer on skymap pixels
457  const double* pixels = m_map.pixels();
458 
459  // Sum event cube
460  if (size() > 0 && pixels != NULL) {
461  for (int i = 0; i < size(); ++i) {
462  number += pixels[i];
463  }
464  }
465 
466  // Return
467  return int(number+0.5);
468 }
469 
470 
471 /***********************************************************************//**
472  * @brief Set event cube from sky map
473  ***************************************************************************/
474 void GLATEventCube::map(const GSkyMap& map)
475 {
476  // Store sky map
477  m_map = map;
478 
479  // Compute sky directions
480  set_directions();
481 
482  // Return
483  return;
484 }
485 
486 
487 /***********************************************************************//**
488  * @brief Print event cube information
489  *
490  * @param[in] chatter Chattiness.
491  * @return String containing event cube information.
492  ***************************************************************************/
493 std::string GLATEventCube::print(const GChatter& chatter) const
494 {
495  // Initialise result string
496  std::string result;
497 
498  // Continue only if chatter is not silent
499  if (chatter != SILENT) {
500 
501  // Append header
502  result.append("=== GLATEventCube ===");
503 
504  // Append information
505  result.append("\n"+gammalib::parformat("Number of elements") +
506  gammalib::str(size()));
507  result.append("\n"+gammalib::parformat("Number of pixels"));
508  result.append(gammalib::str(m_map.nx()) +
509  " x " +
510  gammalib::str(m_map.ny()));
511  result.append("\n"+gammalib::parformat("Number of energy bins") +
512  gammalib::str(ebins()));
513  result.append("\n"+gammalib::parformat("Number of events") +
514  gammalib::str(number()));
515 
516  // Append time interval
517  result.append("\n"+gammalib::parformat("Time interval"));
518  if (gti().size() > 0) {
519  result.append(gammalib::str(tstart().secs()) +
520  " - " +
521  gammalib::str(tstop().secs())+" sec");
522  }
523  else {
524  result.append("not defined");
525  }
526 
527  // Append energy range
528  result.append("\n"+gammalib::parformat("Energy range"));
529  if (ebounds().size() > 0) {
530  result.append(emin().print()+" - "+emax().print(chatter));
531  }
532  else {
533  result.append("not defined");
534  }
535 
536  // Append detailed information
537  GChatter reduced_chatter = gammalib::reduce(chatter);
538  if (reduced_chatter > SILENT) {
539 
540  // Append sky projection
541  if (m_map.projection() != NULL) {
542  result.append("\n"+m_map.projection()->print(reduced_chatter));
543  }
544 
545  // Append source maps
546  result.append("\n"+gammalib::parformat("Number of source maps"));
547  result.append(gammalib::str(m_srcmap.size()));
548  for (int i = 0; i < m_srcmap.size(); ++i) {
549  result.append("\n"+gammalib::parformat(" "+m_srcmap_names[i]));
550  result.append(gammalib::str(m_srcmap[i]->nx()));
551  result.append(" x ");
552  result.append(gammalib::str(m_srcmap[i]->ny()));
553  result.append(" x ");
554  result.append(gammalib::str(m_srcmap[i]->nmaps()));
555  }
556 
557  }
558 
559  } // endif: chatter was not silent
560 
561  // Return result
562  return result;
563 }
564 
565 
566 /***********************************************************************//**
567  * @brief Return name of diffuse model
568  *
569  * @param[in] index Diffuse model index [0,...,ndiffrsp()-1].
570  * @return Name of diffuse model.
571  *
572  * @exception GException::out_of_range
573  * Model index out of valid range.
574  *
575  * Returns name of diffuse model.
576  ***************************************************************************/
577 std::string GLATEventCube::diffname(const int& index) const
578 {
579  // Optionally check if the index is valid
580  #if defined(G_RANGE_CHECK)
581  if (index < 0 || index >= ndiffrsp()) {
582  throw GException::out_of_range(G_DIFFNAME, "Diffuse model index",
583  index, ndiffrsp());
584  }
585  #endif
586 
587  // Return
588  return m_srcmap_names[index];
589 }
590 
591 
592 /***********************************************************************//**
593  * @brief Return diffuse response map
594  *
595  * @param[in] index Diffuse model index [0,...,ndiffrsp()-1].
596  * @return Pointer to diffuse response map.
597  *
598  * @exception GException::out_of_range
599  * Model index out of valid range.
600  *
601  * Returns pointer to diffuse model sky map.
602  ***************************************************************************/
603 GSkyMap* GLATEventCube::diffrsp(const int& index) const
604 {
605  // Optionally check if the index is valid
606  #if defined(G_RANGE_CHECK)
607  if (index < 0 || index >= ndiffrsp()) {
608  throw GException::out_of_range(G_DIFFRSP, "Diffuse model index",
609  index, ndiffrsp());
610  }
611  #endif
612 
613  // Return
614  return m_srcmap[index];
615 }
616 
617 
618 /***********************************************************************//**
619  * @brief Computes the maximum radius (in degrees) around a given source
620  * direction that fits spatially into the event cube
621  *
622  * @param[in] srcDir Source direction.
623  * @return Maximum radius in degrees that fully fits into event cube.
624  *
625  * By computing the sky directions of the event cube boundaries, the maximum
626  * radius is computed that fits fully within the event cube. This method is
627  * used for PSF normalization.
628  ***************************************************************************/
629 double GLATEventCube::maxrad(const GSkyDir& srcDir) const
630 {
631  // Initialise radius
632  double radius = 0.0;
633 
634  // Continue only if sky direction is within sky map
635  if (m_map.contains(srcDir)) {
636 
637  // Set to largest possible radius
638  radius = 180.0;
639 
640  // Move along upper edge in longitude
641  int iy = 0;
642  for (int ix = 0; ix < nx(); ++ix) {
643  GSkyPixel pixel = GSkyPixel(double(ix), double(iy));
644  double distance = m_map.pix2dir(pixel).dist_deg(srcDir);
645  if (distance < radius) {
646  radius = distance;
647  }
648  }
649 
650  // Move along lower edge in longitude
651  iy = ny()-1;
652  for (int ix = 0; ix < nx(); ++ix) {
653  GSkyPixel pixel = GSkyPixel(double(ix), double(iy));
654  double distance = m_map.pix2dir(pixel).dist_deg(srcDir);
655  if (distance < radius) {
656  radius = distance;
657  }
658  }
659 
660  // Move along left edge in latitude
661  int ix = 0;
662  for (int iy = 0; iy < ny(); ++iy) {
663  GSkyPixel pixel = GSkyPixel(double(ix), double(iy));
664  double distance = m_map.pix2dir(pixel).dist_deg(srcDir);
665  if (distance < radius) {
666  radius = distance;
667  }
668  }
669 
670  // Move along right edge in latitude
671  ix = nx()-1;
672  for (int iy = 0; iy < ny(); ++iy) {
673  GSkyPixel pixel = GSkyPixel(double(ix), double(iy));
674  double distance = m_map.pix2dir(pixel).dist_deg(srcDir);
675  if (distance < radius) {
676  radius = distance;
677  }
678  }
679 
680  } // endif: sky direction within sky map
681 
682  // Return radius
683  return radius;
684 }
685 
686 
687 /*==========================================================================
688  = =
689  = Private methods =
690  = =
691  ==========================================================================*/
692 
693 /***********************************************************************//**
694  * @brief Initialise class members
695  ***************************************************************************/
697 {
698  // Initialise members
699  m_bin.clear();
700  m_map.clear();
701  m_time.clear();
702  m_srcmap.clear();
703  m_srcmap_names.clear();
704  m_enodes.clear();
705  m_dirs.clear();
706  m_solidangle.clear();
707  m_energies.clear();
708  m_ewidth.clear();
709  m_ontime = 0.0;
710 
711  // Return
712  return;
713 }
714 
715 
716 /***********************************************************************//**
717  * @brief Copy class members
718  *
719  * @param[in] cube Fermi/LAT event cube.
720  ***************************************************************************/
722 {
723  // Copy Fermi/LAT specific attributes
724  m_bin = cube.m_bin;
725  m_map = cube.m_map;
726  m_time = cube.m_time;
727  m_ontime = cube.m_ontime;
728  m_enodes = cube.m_enodes;
729  m_dirs = cube.m_dirs;
730  m_solidangle = cube.m_solidangle;
731  m_energies = cube.m_energies;
732  m_ewidth = cube.m_ewidth;
733 
734  // Clone source maps and copy source map names
735  m_srcmap.clear();
736  for (int i = 0; i < cube.m_srcmap.size(); ++i) {
737  m_srcmap.push_back(cube.m_srcmap[i]->clone());
738  }
740 
741  // Return
742  return;
743 }
744 
745 
746 /***********************************************************************//**
747  * @brief Delete class members
748  ***************************************************************************/
750 {
751  // Release source maps
752  for (int i = 0; i < m_srcmap.size(); ++i) {
753  if (m_srcmap[i] != NULL) delete m_srcmap[i];
754  m_srcmap[i] = NULL;
755  }
756  m_srcmap.clear();
757  m_srcmap_names.clear();
758 
759  // Return
760  return;
761 }
762 
763 
764 /***********************************************************************//**
765  * @brief Read Fermi/LAT counts map from HDU.
766  *
767  * @param[in] hdu Image HDU.
768  *
769  * This method reads a Fermi/LAT counts map from a FITS image. The counts map
770  * is stored in a GSkyMap object, and a pointer is set up to access the
771  * pixels individually. Recall that skymap pixels are stored in the order
772  * (ix,iy,ebin).
773  ***************************************************************************/
775 {
776  // Load counts map as sky map
777  m_map.read(hdu);
778 
779  // Set sky directions
780  set_directions();
781 
782  // Return
783  return;
784 }
785 
786 
787 /***********************************************************************//**
788  * @brief Read LAT source map from HDU.
789  *
790  * @param[in] hdu Image HDU.
791  *
792  * @exception GException::invalid_argument
793  * Source map in @p hdu is not compatible with event cube.
794  *
795  * This method reads a LAT source map from a FITS image. The source map is
796  * stored in a GSkyMap object and is given in units of counts/pixel/MeV.
797  ***************************************************************************/
799 {
800  // Allocate skymap
801  GSkyMap* map = new GSkyMap;
802 
803  // Read skymap
804  map->read(hdu);
805 
806  // Check that source map sky projection is consistent with counts
807  // map sky projection
808  if (*(m_map.projection()) != *(map->projection())) {
809  std::string msg = "Specified source map projection in HDU \""+
810  hdu.extname()+"\" is incompatible with map projection "
811  "of the event cube. Please specify a compatible "
812  "map projection.";
814  }
815 
816  // Check that source map dimension is consistent with counts map
817  // dimension
818  if (m_map.nx() != map->nx() ||
819  m_map.ny() != map->ny()) {
820  std::string msg = "Number of pixels ("+gammalib::str(map->nx())+","+
821  gammalib::str(map->ny())+") in source map in HDU \""+
822  hdu.extname()+"\" is not compatible with number of "
823  "pixels ("+gammalib::str(m_map.nx())+","+
824  gammalib::str(m_map.ny())+") in the event cube. "
825  "Please specify a source map of compatible size.";
827  }
828 
829  // Check that source map has required number of energy bins
830  if (m_map.nmaps()+1 != map->nmaps()) {
831  std::string msg = "Number of maps ("+gammalib::str(map->nmaps())+")"+
832  " in source map in HDU \""+hdu.extname()+"\" is not "
833  "compatible with number of energy bin boundaries ("+
834  gammalib::str(m_map.nmaps()+1)+") of the event cube. "
835  "Please specify a source map of compatible size.";
837  }
838 
839  // Append source map to list of maps
840  m_srcmap.push_back(map);
841  m_srcmap_names.push_back(hdu.extname());
842 
843  // Return
844  return;
845 }
846 
847 
848 /***********************************************************************//**
849  * @brief Read energy boundaries from HDU.
850  *
851  * @param[in] hdu Energy boundaries table.
852  *
853  * Read the energy boundaries from the HDU.
854  *
855  * @todo Energy bounds read method should take const GFitsTable* as argument
856  ***************************************************************************/
858 {
859  // Read energy boundaries
860  m_ebounds.read(hdu);
861 
862  // Set log mean energies and energy widths
863  set_energies();
864 
865  // Return
866  return;
867 }
868 
869 
870 /***********************************************************************//**
871  * @brief Read GTIs from HDU.
872  *
873  * @param[in] hdu GTI table.
874  *
875  * Reads the Good Time Intervals from the GTI extension. Since the Fermi
876  * LAT Science Tools do not set corrently the time reference for source
877  * maps, the method automatically adds this missing information so that
878  * the time reference is set correctly. The time reference that is assumed
879  * for Fermi LAT is
880  *
881  * MJDREFI 51910
882  * MJDREFF 0.00074287037037037
883  *
884  ***************************************************************************/
886 {
887  // Work on a local copy of the HDU to make the kluge work
888  GFitsTable* hdu_local = hdu.clone();
889 
890  // Kluge: modify HDU table in case that the MJDREF header keyword is
891  // blank. This happens for Fermi LAT source maps since the Science
892  // Tools do not properly write out the time reference. We hard-code
893  // here the Fermi LAT time reference to circumvent the problem.
894  if (hdu.has_card("MJDREF")) {
895  if (gammalib::strip_whitespace(hdu.string("MJDREF")).empty()) {
896  const_cast<GFitsHeader&>(hdu_local->header()).remove("MJDREF");
897  hdu_local->card("MJDREFI", 51910,
898  "Integer part of MJD reference");
899  hdu_local->card("MJDREFF", 0.00074287037037037,
900  "Fractional part of MJD reference");
901  }
902  }
903 
904  // Read Good Time Intervals
905  m_gti.read(*hdu_local);
906 
907  // Set time
908  set_times();
909 
910  // Return
911  return;
912 }
913 
914 
915 /***********************************************************************//**
916  * @brief Set sky directions and solid angles of events cube.
917  *
918  * @exception GException::invalid_value
919  * No sky pixels found in event cube.
920  *
921  * This method computes the sky directions and solid angles for all event
922  * cube pixels. Sky directions are stored in an array of GLATInstDir objects
923  * while solid angles are stored in units of sr in a double precision array.
924  ***************************************************************************/
926 {
927  // Throw an exception if we have no sky pixels
928  if (npix() < 1) {
929  std::string msg = "LAT event cube contains no sky pixels. Please "
930  "provide a valid LAT event cube.";
932  }
933 
934  // Clear old pixel directions and solid angle
935  m_dirs.clear();
936  m_solidangle.clear();
937 
938  // Reserve space for pixel directions and solid angles
939  m_dirs.reserve(npix());
940  m_solidangle.reserve(npix());
941 
942  // Set pixel directions and solid angles
943  for (int iy = 0; iy < ny(); ++iy) {
944  for (int ix = 0; ix < nx(); ++ix) {
945  GSkyPixel pixel = GSkyPixel(double(ix), double(iy));
946  m_dirs.push_back(GLATInstDir(m_map.pix2dir(pixel)));
947  m_solidangle.push_back(m_map.solidangle(pixel));
948  }
949  }
950 
951  // Return
952  return;
953 }
954 
955 
956 /***********************************************************************//**
957  * @brief Set log mean energies and energy widths of event cube.
958  *
959  * @exception GException::invalid_value
960  * No energy boundaries found in event cube.
961  *
962  * This method computes the log mean energies and the energy widths of the
963  * event cube. The log mean energies and energy widths are stored unit
964  * independent in arrays of GEnergy objects.
965  ***************************************************************************/
967 {
968  // Throw an error if we have no energy bins
969  if (ebins() < 1) {
970  std::string msg = "LAT event cube contains no energy bins. Please "
971  "provide a valid LAT event cube.";
973  }
974 
975  // Clear old bin energies and energy widths
976  m_energies.clear();
977  m_ewidth.clear();
978  m_enodes.clear();
979 
980  // Reserve space for bin energies and energy widths
981  m_energies.reserve(ebins());
982  m_ewidth.reserve(ebins());
983 
984  // Setup bin energies, energy widths and energy nodes
985  for (int i = 0; i < ebins(); ++i) {
986  m_energies.push_back(ebounds().elogmean(i));
987  m_ewidth.push_back(ebounds().emax(i) - ebounds().emin(i));
988  m_enodes.append(log10(ebounds().emin(i).MeV()));
989  }
990  m_enodes.append(log10(ebounds().emax(ebins()-1).MeV()));
991 
992  // Return
993  return;
994 }
995 
996 
997 /***********************************************************************//**
998  * @brief Set mean event time and ontime of event cube
999  *
1000  * @exception GException::invalid_value
1001  * No Good Time Intervals found.
1002  *
1003  * Computes the mean time of the event cube by taking the mean between start
1004  * and stop time. Computes also the ontime by summing up of all good time
1005  * intervals.
1006  *
1007  * @todo Could add a more sophisticated mean event time computation that
1008  * weights by the length of the GTIs, yet so far we do not really use
1009  * the mean event time, hence there is no rush to implement this.
1010  ***************************************************************************/
1012 {
1013  // Throw an error if GTI is empty
1014  if (m_gti.size() < 1) {
1015  std::string msg = "No Good Time Intervals have been found in event "
1016  "cube. Every LAT event cube needs a definition "
1017  "of the Good Time Intervals.";
1019  }
1020 
1021  // Compute mean time
1022  m_time = m_gti.tstart() + 0.5 * (m_gti.tstop() - m_gti.tstart());
1023 
1024  // Set ontime
1025  m_ontime = m_gti.ontime();
1026 
1027  // Return
1028  return;
1029 }
1030 
1031 
1032 /***********************************************************************//**
1033  * @brief Set event bin
1034  *
1035  * @param[in] index Event index [0,...,size()-1].
1036  *
1037  * @exception GException::out_of_range
1038  * Event index is outside valid range.
1039  * @exception GException::invalid_value
1040  * Energy vectors have not been set up.
1041  * Sky directions and solid angles vectors have not been set up.
1042  *
1043  * This method provides the event attributes to the event bin. The event bin
1044  * is in fact physically stored in the event cube, and only a single event
1045  * bin is indeed allocated. This method sets up the pointers in the event
1046  * bin so that a client can easily access the information of individual bins
1047  * as if they were stored in an array.
1048  ***************************************************************************/
1049 void GLATEventCube::set_bin(const int& index)
1050 {
1051  // Optionally check if the index is valid
1052  #if defined(G_RANGE_CHECK)
1053  if (index < 0 || index >= size()) {
1054  throw GException::out_of_range(G_SET_BIN, "Event index", index, size());
1055  }
1056  #endif
1057 
1058  // Check for the existence of energies and energy widths
1059  if (m_energies.size() != ebins() || m_ewidth.size() != ebins()) {
1060  std::string msg = "Number of energy bins ("+gammalib::str(ebins())+
1061  ") is incompatible with size of energies ("+
1062  gammalib::str(m_energies.size())+") an/or size of "
1063  "energy widths ("+gammalib::str(m_ewidth.size())+
1064  "). Please provide an correctly initialised event "
1065  "cube.";
1067  }
1068 
1069  // Check for the existence of sky directions and solid angles
1070  if (m_dirs.size() != npix() || m_solidangle.size() != npix()) {
1071  std::string msg = "Number of sky pixels ("+gammalib::str(npix())+
1072  ") is incompatible with size of sky directions ("+
1073  gammalib::str(m_dirs.size())+") and/or size of "
1074  "solid angles ("+gammalib::str(m_solidangle.size())+
1075  "). Please provide an correctly initialised event "
1076  "cube.";
1078  }
1079 
1080  // Get pixel and energy bin indices.
1081  m_bin.m_index = index;
1082  m_bin.m_ipix = index % npix();
1083  m_bin.m_ieng = index / npix();
1084 
1085  // Set pointers
1086  m_bin.m_cube = this;
1087  m_bin.m_counts = const_cast<double*>(&(m_map.pixels()[index]));
1089  m_bin.m_time = &m_time;
1090  m_bin.m_dir = &(m_dirs[m_bin.m_ipix]);
1093  m_bin.m_ontime = &m_ontime;
1094 
1095  // Return
1096  return;
1097 }
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:286
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:482
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
const GSkyProjection * projection(void) const
Returns pointer to sky projection.
Definition: GSkyMap.hpp:463
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:48
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:252
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:761
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:80
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
int size(void) const
Return number of Good Time Intervals.
Definition: GGti.hpp:154
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:1293
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:1360
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:389
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.
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
#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:2023
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:1100
Fermi/LAT event cube class.
const GTime & tstop(void) const
Returns latest stop time in Good Time Intervals.
Definition: GGti.hpp:207
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:194
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:368
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:753
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.
Exception handler interface definition.
#define G_DIFFNAME
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
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:1858
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:1143
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:345
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: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
int ny(void) const
Return number of bins in Y direction.
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
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:1295
GEnergy elogmean(const GEnergy &a, const GEnergy &b)
Computes log mean energy.
Definition: GTools.cpp:1290
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:489
FITS table abstract base class interface definition.