GammaLib  1.7.0.dev
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
GCTACubeEdisp.cpp
Go to the documentation of this file.
1 /***************************************************************************
2  * GCTACubeEdisp.cpp - CTA cube analysis energy dispersion class *
3  * ----------------------------------------------------------------------- *
4  * copyright (C) 2016-2018 by Michael Mayer *
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 GCTACubeEdisp.cpp
23  * @brief CTA cube analysis energy dispersion class implementation
24  * @author Michael Mayer
25  */
26 
27 /* __ Includes ___________________________________________________________ */
28 #ifdef HAVE_CONFIG_H
29 #include <config.h>
30 #endif
31 #include "GTools.hpp"
32 #include "GMath.hpp"
33 #include "GLog.hpp"
34 #include "GFits.hpp"
35 #include "GFitsImage.hpp"
36 #include "GFitsTable.hpp"
37 #include "GObservations.hpp"
38 #include "GSkyRegionCircle.hpp"
39 #include "GCTACubeEdisp.hpp"
40 #include "GCTAObservation.hpp"
41 #include "GCTAResponseIrf.hpp"
42 #include "GCTAEventList.hpp"
43 #include "GCTAEventCube.hpp"
44 
45 /* __ Method name definitions ____________________________________________ */
46 #define G_CONSTRUCTOR1 "GCTACubeEdisp(GCTAEventCube&, double&, int&)"
47 #define G_CONSTRUCTOR2 "GCTACubeEdisp(std::string&, std::string&, double&, "\
48  "double&, double&, double&, int&, int&, GEbounds&, "\
49  "double&, int&)"
50 #define G_EBOUNDS "GCTACubeEdisp::ebounds(GEnergy&)"
51 #define G_FILL_CUBE "GCTACubeEdisp::fill_cube(GCTAObservation&, GSkyMap*, "\
52  "GLog*)"
53 
54 /* __ Macros _____________________________________________________________ */
55 
56 /* __ Coding definitions _________________________________________________ */
57 
58 /* __ Debug definitions __________________________________________________ */
59 
60 /* __ Constants __________________________________________________________ */
61 
62 
63 /*==========================================================================
64  = =
65  = Constructors/destructors =
66  = =
67  ==========================================================================*/
68 
69 /***********************************************************************//**
70  * @brief Void constructor
71  *
72  * Constructs an empty energy dispersion cube.
73  ***************************************************************************/
75 {
76  // Initialise class members
77  init_members();
78 
79  // Return
80  return;
81 }
82 
83 
84 /***********************************************************************//**
85  * @brief Copy constructor
86  *
87  * @param[in] cube Energy dispersion cube.
88  *
89  * Constructs an energy dispersion cube by copying another energy dispersion
90  * cube.
91  ***************************************************************************/
93 {
94  // Initialise class members
95  init_members();
96 
97  // Copy members
98  copy_members(cube);
99 
100  // Return
101  return;
102 }
103 
104 
105 /***********************************************************************//**
106  * @brief File constructor
107  *
108  * @param[in] filename Energy dispersion cube filename.
109  *
110  * Constructs an energy dispersion cube by loading the energy dispersion
111  * information from an energy dispersion cube FITS file. See the load()
112  * method for details about the format of the FITS file.
113  ***************************************************************************/
115 {
116  // Initialise class members
117  init_members();
118 
119  // Load Edisp cube from file
120  load(filename);
121 
122  // Return
123  return;
124 }
125 
126 
127 /***********************************************************************//**
128  * @brief Event cube constructor
129  *
130  * @param[in] cube Event cube.
131  * @param[in] mmax Maximum energy migration.
132  * @param[in] nmbins Number of migration bins (2 ...).
133  *
134  * @exception GException::invalid_argument
135  * Maximum energy migration or number of migration bins invalid.
136  *
137  * Construct an energy dispersion cube with all elements set to zero using
138  * the same binning and sky projection that is used for an event cube.
139  ***************************************************************************/
140 GCTACubeEdisp::GCTACubeEdisp(const GCTAEventCube& cube, const double& mmax,
141  const int& nmbins)
142 {
143  // Throw an exception if the number of migra bins is invalid or the
144  // maximum migration is not positive
145  if (nmbins < 2) {
146  std::string msg = "Number "+gammalib::str(nmbins)+" of migration "
147  "bins is smaller than 2. Please request at least "
148  "2 migration bins.";
150  }
151  if (mmax <= 0.0) {
152  std::string msg = "Maximum migration "+gammalib::str(mmax)+" is not "
153  "positive. Please specify a positive maximum "
154  "energy migration.";
156  }
157 
158  // Initialise class members
159  init_members();
160 
161  // Store energy boundaries
162  m_energies.set(cube.ebounds());
163 
164  // Set energy node array used for interpolation
165  set_eng_axis();
166 
167  // Set migration axis used for interpolation
168  set_migras(mmax, nmbins);
169 
170  // Compute number of sky maps
171  int nmaps = m_energies.size() * m_migras.size();
172 
173  // Set energy dispersion cube to event cube
174  m_cube = cube.counts();
175 
176  // Set appropriate number of skymaps
177  m_cube.nmaps(nmaps);
178 
179  // Set cube shape
181 
182  // Set all energy dispersion cube pixels to zero as we want to have
183  // a clean map upon construction
184  m_cube = 0.0;
185 
186  // Return
187  return;
188 }
189 
190 
191 /***********************************************************************//**
192  * @brief Energy dispersion cube constructor
193  *
194  * @param[in] wcs World Coordinate System.
195  * @param[in] coords Coordinate System (CEL or GAL).
196  * @param[in] x X coordinate of sky map centre (deg).
197  * @param[in] y Y coordinate of sky map centre (deg).
198  * @param[in] dx Pixel size in x direction at centre (deg/pixel).
199  * @param[in] dy Pixel size in y direction at centre (deg/pixel).
200  * @param[in] nx Number of pixels in x direction.
201  * @param[in] ny Number of pixels in y direction.
202  * @param[in] energies True energies.
203  * @param[in] mmax Maximum energy migration.
204  * @param[in] nmbins Number of migration bins (2 ...).
205  *
206  * @exception GException::invalid_argument
207  * Maximum energy migration or number of migration bins invalid.
208  *
209  * Constructs an energy dispersion cube by specifying the sky map grid and
210  * the energies.
211  ***************************************************************************/
212 GCTACubeEdisp::GCTACubeEdisp(const std::string& wcs,
213  const std::string& coords,
214  const double& x,
215  const double& y,
216  const double& dx,
217  const double& dy,
218  const int& nx,
219  const int& ny,
220  const GEnergies& energies,
221  const double& mmax,
222  const int& nmbins)
223 {
224  // Throw an exception if the number of migra bins is invalid or the
225  // maximum migration is not positive
226  if (nmbins < 2) {
227  std::string msg = "Number "+gammalib::str(nmbins)+" of migration "
228  "bins is smaller than 2. Please request at least "
229  "2 migration bins.";
231  }
232  if (mmax <= 0.0) {
233  std::string msg = "Maximum migration "+gammalib::str(mmax)+" is not "
234  "positive. Please specify a positive maximum "
235  "energy migration.";
237  }
238 
239  // Initialise class members
240  init_members();
241 
242  // Store true energies
244 
245  // Set energy node array used for interpolation
246  set_eng_axis();
247 
248  // Set migration axis used for interpolation
249  set_migras(mmax, nmbins);
250 
251  // Compute number of sky maps
252  int nmaps = m_energies.size() * m_migras.size();
253 
254  // Create sky map
255  m_cube = GSkyMap(wcs, coords, x, y, dx, dy, nx, ny, nmaps);
256 
257  // Set cube shape
259 
260  // Return
261  return;
262 }
263 
264 
265 /***********************************************************************//**
266  * @brief Destructor
267  *
268  * Destructs energy dispersion cube.
269  ***************************************************************************/
271 {
272  // Free members
273  free_members();
274 
275  // Return
276  return;
277 }
278 
279 
280 /*==========================================================================
281  = =
282  = Operators =
283  = =
284  ==========================================================================*/
285 
286 /***********************************************************************//**
287  * @brief Assignment operator
288  *
289  * @param[in] cube Energy dispersion cube.
290  * @return Energy dispersion cube.
291  *
292  * Assigns energy dispersion cube.
293  ***************************************************************************/
295 {
296  // Execute only if object is not identical
297  if (this != &cube) {
298 
299  // Free members
300  free_members();
301 
302  // Initialise private members
303  init_members();
304 
305  // Copy members
306  copy_members(cube);
307 
308  } // endif: object was not identical
309 
310  // Return this object
311  return *this;
312 }
313 
314 
315 /***********************************************************************//**
316  * @brief Return energy dispersion in units of MeV\f$^{-1}\f$
317  *
318  * @param[in] ereco Reconstructed event energy.
319  * @param[in] etrue True photon energy.
320  * @param[in] dir Coordinate of the true photon position.
321  * @return Energy dispersion (MeV\f$^{-1}\f$)
322  *
323  * Returns the energy dispersion for a given energy migration (reconstructed
324  * over true energy), true photon energy, and sky direction in units of
325  * MeV\f$^{-1}\f$.
326  ***************************************************************************/
327 double GCTACubeEdisp::operator()(const GEnergy& ereco,
328  const GEnergy& etrue,
329  const GSkyDir& dir) const
330 {
331  // Update indices and weighting factors for interpolation
332  update(ereco, etrue);
333 
334  // Perform bi-linear interpolation
335  double edisp = m_wgt1 * m_cube(dir, m_inx1) +
336  m_wgt2 * m_cube(dir, m_inx2) +
337  m_wgt3 * m_cube(dir, m_inx3) +
338  m_wgt4 * m_cube(dir, m_inx4);
339 
340  // Make sure that energy dispersion does not become negative
341  if (edisp < 0.0) {
342  edisp = 0.0;
343  }
344 
345  // Return energy dispersion
346  return edisp;
347 }
348 
349 
350 /*==========================================================================
351  = =
352  = Public methods =
353  = =
354  ==========================================================================*/
355 
356 /***********************************************************************//**
357  * @brief Clear energy dispersion cube
358  *
359  * Clears energy dispersion cube.
360  ***************************************************************************/
362 {
363  // Free class members
364  free_members();
365 
366  // Initialise members
367  init_members();
368 
369  // Return
370  return;
371 }
372 
373 
374 /***********************************************************************//**
375  * @brief Clone energy dispersion cube
376  *
377  * @return Pointer to deep copy of energy dispersion cube.
378  ***************************************************************************/
380 {
381  return new GCTACubeEdisp(*this);
382 }
383 
384 
385 /***********************************************************************//**
386  * @brief Set energy dispersion cube for one CTA observation
387  *
388  * @param[in] obs CTA observation.
389  *
390  * Sets the energy dispersion cube for one CTA observation.
391  ***************************************************************************/
393 {
394  // Clear energy dispersion cube
395  clear_cube();
396 
397  // Fill energy dispersion cube
398  fill_cube(obs);
399 
400  // Return
401  return;
402 }
403 
404 
405 /***********************************************************************//**
406  * @brief Fill energy dispersion cube from observation container
407  *
408  * @param[in] obs Observation container.
409  * @param[in] log Pointer towards logger.
410  *
411  * Sets the energy dispersion cube from all CTA observations in the
412  * observation container.
413  ***************************************************************************/
415 {
416  // Clear energy dispersion cube
417  clear_cube();
418 
419  // Initialise skymap for exposure weight accumulation
420  GSkyMap exposure(m_cube);
421 
422  // Loop over all observations in container
423  for (int i = 0; i < obs.size(); ++i) {
424 
425  // Get observation and continue only if it is a CTA observation
426  const GCTAObservation* cta = dynamic_cast<const GCTAObservation*>
427  (obs[i]);
428 
429  // Skip observation if it's not CTA
430  if (cta == NULL) {
431  if (log != NULL) {
432  *log << "Skipping ";
433  *log << obs[i]->instrument();
434  *log << " observation ";
435  *log << "\"" << obs[i]->name() << "\"";
436  *log << " (id=" << obs[i]->id() << ")" << std::endl;
437  }
438  continue;
439  }
440 
441  // Skip observation if we have a binned observation
442  if (cta->eventtype() == "CountsCube") {
443  if (log != NULL) {
444  *log << "Skipping binned ";
445  *log << cta->instrument();
446  *log << " observation ";
447  *log << "\"" << cta->name() << "\"";
448  *log << " (id=" << cta->id() << ")" << std::endl;
449  }
450  continue;
451  }
452 
453  // Fill exposure cube cube
454  fill_cube(*cta, &exposure, log);
455 
456  } // endfor: looped over observations
457 
458  // Compute mean energy dispersion cube by dividing through the weights
459  for (int pixel = 0; pixel < m_cube.npix(); ++pixel) {
460  for (int iebin = 0; iebin < m_energies.size(); ++iebin) {
461  if (exposure(pixel, iebin) > 0.0) {
462  double norm = 1.0 / exposure(pixel, iebin);
463  for (int imigra = 0; imigra < m_migras.size(); ++imigra) {
464  int imap = offset(imigra, iebin);
465  m_cube(pixel, imap) *= norm;
466  }
467  }
468  else {
469  for (int imigra = 0; imigra < m_migras.size(); ++imigra) {
470  int imap = offset(imigra, iebin);
471  m_cube(pixel, imap) = 0.0;
472  }
473  }
474  }
475  }
476 
477  // Return
478  return;
479 }
480 
481 
482 /***********************************************************************//**
483  * @brief Read energy dispersion cube from FITS file
484  *
485  * @param[in] fits FITS file.
486  *
487  * Reads the energy dispersion cube from a FITS file. The energy dispersion
488  * cube values are expected in the primary extension of the FITS file, the
489  * true energy boundaries are expected in the `ENERGIES` extension, and the
490  * migration values are expected in the `MIGRAS` extension.
491  ***************************************************************************/
492 void GCTACubeEdisp::read(const GFits& fits)
493 {
494  // Clear object
495  clear();
496 
497  // Get HDUs
498  const GFitsImage& hdu_edispcube = *fits.image("Primary");
499  const GFitsTable& hdu_energies = *fits.table(gammalib::extname_energies);
500  const GFitsTable& hdu_migras = *fits.table(gammalib::extname_cta_migras);
501 
502  // Read energy dispersion cube
503  m_cube.read(hdu_edispcube);
504 
505  // Read true energy nodes
506  m_energies.read(hdu_energies);
507 
508  // Read migration nodes
509  m_migras.read(hdu_migras);
510 
511  // Set energy node array used for interpolation
512  set_eng_axis();
513 
514  // Compute true energy boundaries
515  compute_ebounds();
516 
517  // Return
518  return;
519 }
520 
521 
522 /***********************************************************************//**
523  * @brief Write energy dispersion cube into FITS file.
524  *
525  * @param[in] fits FITS file.
526  *
527  * Write the energy dispersion cube into a FITS file. The energy dispersion
528  * cube values are written into the primary extension of the FITS file, the
529  * true energy boundaries are written into the `ENERGIES` extension, and the
530  * migration values are written into the `MIGRAS` extension.
531  ***************************************************************************/
532 void GCTACubeEdisp::write(GFits& fits) const
533 {
534  // Remove extensions from FITS file
537  }
540  }
541  if (fits.contains("Primary")) {
542  fits.remove("Primary");
543  }
544 
545  // Write cube
546  m_cube.write(fits);
547 
548  // Write energy boundaries
550 
551  // Write migration nodes
553 
554  // Return
555  return;
556 }
557 
558 
559 /***********************************************************************//**
560  * @brief Load energy dispersion cube from FITS file
561  *
562  * @param[in] filename FITS file name.
563  *
564  * Loads the energy dispersion cube from a FITS file. See the read() method
565  * for information about the expected FITS file structure.
566  ***************************************************************************/
567 void GCTACubeEdisp::load(const GFilename& filename)
568 {
569  // Put into OpenMP criticial zone
570  #pragma omp critical(GCTACubeEdisp_load)
571  {
572 
573  // Open FITS file
574  GFits fits(filename);
575 
576  // Read Edisp cube
577  read(fits);
578 
579  // Close Edisp file
580  fits.close();
581 
582  } // end of OpenMP critical zone
583 
584  // Store filename
586 
587  // Return
588  return;
589 }
590 
591 
592 /***********************************************************************//**
593  * @brief Save energy dispersion cube into FITS file
594  *
595  * @param[in] filename FITS file name.
596  * @param[in] clobber Overwrite existing file?
597  *
598  * Save the energy dispersion cube into a FITS file.
599  ***************************************************************************/
600 void GCTACubeEdisp::save(const GFilename& filename, const bool& clobber) const
601 {
602  // Put into OpenMP criticial zone
603  #pragma omp critical(GCTACubeEdisp_save)
604  {
605 
606  // Create FITS file
607  GFits fits;
608 
609  // Write PSF cube
610  write(fits);
611 
612  // Save FITS file
613  fits.saveto(filename, clobber);
614 
615  // Close Edisp file
616  fits.close();
617 
618  } // end of OpenMP critical zone
619 
620  // Store filename
622 
623  // Return
624  return;
625 }
626 
627 
628 /***********************************************************************//**
629  * @brief Return boundaries in true energy
630  *
631  * @param[in] obsEng Observed photon energy.
632  * @return Boundaries in true energy.
633  *
634  * @exception GException::invalid_value
635  * No energies defined for energy dispersion cube
636  *
637  * Returns the boundaries in true photon energies that enclose the energy
638  * disperson for a given observed photon energy @p obsEng.
639  ***************************************************************************/
641 {
642  // Throw an exception if there are no energies
643  if (m_energies.is_empty()) {
644  std::string msg = "No energies defined for energy dispersion cube. "
645  "Please define the energies before calling the "
646  "method.";
648  }
649 
650  // If ebounds were not computed, before, compute them now
651  if (m_ebounds.empty()) {
652  compute_ebounds();
653  }
654 
655  // Return the index of the energy boundaries that are just below the
656  // observed energy. As the energy dispersion decreases with increasing
657  // energy we assure in that way that we integrate over a sufficiently
658  // large true energy range.
659  int index = 0;
660  if (obsEng > m_energies[0]) {
661  for (int i = 1; i < m_energies.size(); ++i) {
662  if (obsEng < m_energies[i]) {
663  break;
664  }
665  index++;
666  }
667  if (index >= m_energies.size()) {
668  index = m_energies.size();
669  }
670  }
671 
672  // Return true energy boundaries
673  return (m_ebounds[index]);
674 }
675 
676 
677 /***********************************************************************//**
678  * @brief Print energy dispersion cube information
679  *
680  * @param[in] chatter Chattiness.
681  * @return String containing energy dispersion cube information.
682  ***************************************************************************/
683 std::string GCTACubeEdisp::print(const GChatter& chatter) const
684 {
685  // Initialise result string
686  std::string result;
687 
688  // Continue only if chatter is not silent
689  if (chatter != SILENT) {
690 
691  // Append header
692  result.append("=== GCTACubeEdisp ===");
693 
694  // Append information
695  result.append("\n"+gammalib::parformat("Filename")+m_filename);
696 
697  // Append energies
698  if (m_energies.size() > 0) {
699  result.append("\n"+m_energies.print(chatter));
700  }
701  else {
702  result.append("\n"+gammalib::parformat("Energies") +
703  "Not defined");
704  }
705 
706  // Append number of migra bins
707  result.append("\n"+gammalib::parformat("Number of migra bins") +
709 
710  // Append migra range
711  result.append("\n"+gammalib::parformat("Migra range"));
712  if (m_migras.size() > 0) {
713  result.append(gammalib::str(m_migras[0]));
714  result.append(" - ");
715  result.append(gammalib::str(m_migras[m_migras.size()-1]));
716  }
717  else {
718  result.append("not defined");
719  }
720 
721  // Append skymap definition
722  result.append("\n"+m_cube.print(chatter));
723 
724  } // endif: chatter was not silent
725 
726  // Return result
727  return result;
728 }
729 
730 
731 /*==========================================================================
732  = =
733  = Private methods =
734  = =
735  ==========================================================================*/
736 
737 /***********************************************************************//**
738  * @brief Initialise class members
739  ***************************************************************************/
741 {
742  // Initialise members
743  m_filename.clear();
744  m_cube.clear();
745  m_energies.clear();
746  m_elogmeans.clear();
747  m_migras.clear();
748 
749  // Initialise cache
750  m_inx1 = 0;
751  m_inx2 = 0;
752  m_inx3 = 0;
753  m_inx4 = 0;
754  m_wgt1 = 0.0;
755  m_wgt2 = 0.0;
756  m_wgt3 = 0.0;
757  m_wgt4 = 0.0;
758  m_ebounds.clear();
759 
760  // Return
761  return;
762 }
763 
764 
765 /***********************************************************************//**
766  * @brief Copy class members
767  *
768  * @param[in] cube Energy dispersion cube.
769  ***************************************************************************/
771 {
772  // Copy members
773  m_filename = cube.m_filename;
774  m_cube = cube.m_cube;
775  m_energies = cube.m_energies;
776  m_elogmeans = cube.m_elogmeans;
777  m_migras = cube.m_migras;
778 
779  // Copy cache
780  m_inx1 = cube.m_inx1;
781  m_inx2 = cube.m_inx2;
782  m_inx3 = cube.m_inx3;
783  m_inx4 = cube.m_inx4;
784  m_wgt1 = cube.m_wgt1;
785  m_wgt2 = cube.m_wgt2;
786  m_wgt3 = cube.m_wgt3;
787  m_wgt4 = cube.m_wgt4;
788  m_ebounds = cube.m_ebounds;
789 
790  // Return
791  return;
792 }
793 
794 /***********************************************************************//**
795  * @brief Delete class members
796  ***************************************************************************/
798 {
799  // Return
800  return;
801 }
802 
803 
804 /***********************************************************************//**
805  * @brief Clear all pixels in the energy dispersion cube
806  ***************************************************************************/
808 {
809  // Loop over all maps
810  for (int imap = 0; imap < m_cube.nmaps(); ++imap) {
811 
812  // Loop over all pixels in sky map
813  for (int pixel = 0; pixel < m_cube.npix(); ++pixel) {
814 
815  // Reset cube value to zero
816  m_cube(pixel, imap) = 0.0;
817 
818  } // endfor: looped over all pixels
819 
820  } // endfor: looped over maps
821 
822  // Return
823  return;
824 }
825 
826 
827 /***********************************************************************//**
828  * @brief Fill energy dispersion cube from observation container
829  *
830  * @param[in] obs CTA observation.
831  * @param[in] exposure Pointer towards exposure map.
832  * @param[in] log Pointer towards logger.
833  *
834  * @exception GException::invalid_value
835  * No RoI or response found in CTA observation.
836  *
837  * Fills the energy dispersion cube from the observation container.
838  *
839  * If the @p exposure argument is not NULL, the energy dispersion is weighted
840  * by the exposure, and the exposure weighting is added to the exposure map.
841  *
842  * If the @p log argument is not NULL, the method puts information about
843  * exclusion and inclusion of a CTA observation into the energy dispersion
844  * cube.
845  ***************************************************************************/
847  GSkyMap* exposure,
848  GLog* log)
849 {
850  // Only continue if we have an event list
851  if (obs.eventtype() == "EventList") {
852 
853  // Extract pointing direction, energy boundaries and ROI from
854  // observation
855  GSkyDir pnt = obs.pointing().dir();
856  GEbounds obs_ebounds = obs.ebounds();
857  GCTARoi roi = obs.roi();
858 
859  // Check for RoI sanity
860  if (!roi.is_valid()) {
861  std::string msg = "No RoI information found in "+obs.instrument()+
862  " observation \""+obs.name()+"\". Run ctselect "
863  "to specify an RoI for this observation.";
865  }
866 
867  // Convert RoI into a circular region for overlap checking
868  GSkyRegionCircle roi_reg(roi.centre().dir(), roi.radius());
869 
870  // Extract response from observation
871  const GCTAResponseIrf* rsp = dynamic_cast<const GCTAResponseIrf*>
872  (obs.response());
873  if (rsp == NULL) {
874  std::string msg = "No valid instrument response function found in "+
875  obs.instrument()+" observation \""+obs.name()+
876  "\". Please specify the instrument response "
877  "function for this observation.";
879  }
880 
881  // Get energy dispersion component
882  const GCTAEdisp* edisp = rsp->edisp();
883  if (edisp == NULL) {
884  std::string msg = "No energy dispersion rcomponent found in the "
885  "the instrument response of "+
886  obs.instrument()+" observation \""+obs.name()+
887  "\". Please provide an instrument response that "
888  "comprises an energy dispersion component.";
890  }
891 
892  // Skip observation if livetime is zero
893  if (obs.livetime() == 0.0) {
894  if (log != NULL) {
895  *log << "Skipping unbinned ";
896  *log << obs.instrument();
897  *log << " observation ";
898  *log << "\"" << obs.name() << "\"";
899  *log << " (id=" << obs.id() << ") due to zero livetime";
900  *log << std::endl;
901  }
902  }
903 
904  // Skip observation if observation is outside the bounds of the cube
905  else if (!m_cube.overlaps(roi_reg)) {
906  if (log != NULL) {
907  *log << "Skipping unbinned ";
908  *log << obs.instrument();
909  *log << " observation ";
910  *log << "\"" << obs.name() << "\"";
911  *log << " (id=" << obs.id() << ") since it does not overlap ";
912  *log << "with the energy dispersion cube.";
913  *log << std::endl;
914  }
915  }
916 
917  // ... otherwise continue
918  else {
919 
920  // Announce observation usage
921  if (log != NULL) {
922  *log << "Including ";
923  *log << obs.instrument();
924  *log << " observation \"" << obs.name();
925  *log << "\" (id=" << obs.id() << ")";
926  *log << " in energy dispersion cube computation." << std::endl;
927  }
928 
929  // Loop over all pixels in sky map
930  for (int pixel = 0; pixel < m_cube.npix(); ++pixel) {
931 
932  // Get pixel sky direction
933  GSkyDir dir = m_cube.inx2dir(pixel);
934 
935  // Skip pixel if it is outside the RoI
936  if (roi.centre().dir().dist_deg(dir) > roi.radius()) {
937  continue;
938  }
939 
940  // Compute theta angle with respect to pointing direction in
941  // radians
942  double theta = pnt.dist(dir);
943 
944  // Loop over all energy bins
945  for (int iebin = 0; iebin < m_energies.size(); ++iebin) {
946 
947  // Get log10 of true energy in TeV
948  double logEsrc = m_energies[iebin].log10TeV();
949 
950  // Compute exposure weight. If no exposure map is
951  // specified, the weight is one. Otherwise, the weight
952  // is the effective area for this offset angle and
953  // source energy times the livetime of the observation.
954  double weight = 1.0;
955  if (exposure != NULL) {
956  weight = rsp->aeff(theta, 0.0, 0.0, 0.0, logEsrc) *
957  obs.livetime();
958  (*exposure)(pixel, iebin) += weight;
959  }
960 
961  // Loop over migration values
962  for (int imigra = 0; imigra < m_migras.size(); ++imigra) {
963 
964  // Compute energy migration fraction
965  double migra = m_migras[imigra];
966 
967  // Skip migra bin if zero
968  if (migra <= 0.0) {
969  continue;
970  }
971 
972  // Compute reconstructed energy
973  GEnergy Eobs = migra * m_energies[iebin];
974 
975  // Set map index
976  int imap = offset(imigra, iebin);
977 
978  // Add energy dispersion cube value
979  m_cube(pixel, imap) += rsp->edisp(Eobs,
980  m_energies[iebin],
981  theta, 0.0,
982  0.0, 0.0) * weight;
983 
984  } // endfor: looped over migration bins
985 
986  } // endfor: looped over energy bins
987 
988  } // endfor: looped over all pixels
989 
990  } // endelse: livetime was not zero
991 
992  } // endif: observation contained an event list
993 
994  // Return
995  return;
996 }
997 
998 
999 /***********************************************************************//**
1000  * @brief Update energy dispersion cube indices and weights cache
1001  *
1002  * @param[in] ereco Reconstructed event energy.
1003  * @param[in] etrue True photon energy.
1004  *
1005  * Updates the energy dispersion cube indices, stored in the members m_inx1,
1006  * m_inx2, m_inx3, and m_inx4, and weights cache, stored in m_wgt1, m_wgt2,
1007  * m_wgt3, and m_wgt4.
1008  ***************************************************************************/
1009 void GCTACubeEdisp::update(const GEnergy& ereco, const GEnergy& etrue) const
1010 {
1011  // Set node array for energy interpolation
1012  m_elogmeans.set_value(etrue.log10TeV());
1013 
1014  // Compute migration value
1015  double migra = ereco/etrue;
1016 
1017  // Set node array for migra interpolation
1018  m_migras.set_value(migra);
1019 
1020  // Set indices for bi-linear interpolation
1025 
1026  // Set weighting factors for bi-linear interpolation
1031 
1032  // Return
1033  return;
1034 }
1035 
1036 
1037 /***********************************************************************//**
1038  * @brief Set nodes for interpolation in true energy
1039  ***************************************************************************/
1041 {
1042  // Clear nodes
1043  m_elogmeans.clear();
1044 
1045  // Set nodes
1046  for (int i = 0; i < m_energies.size(); ++i) {
1047  m_elogmeans.append(m_energies[i].log10TeV());
1048  }
1049 
1050  // Return
1051  return;
1052 }
1053 
1054 
1055 /***********************************************************************//**
1056  * @brief Set nodes for interpolation in migration
1057  *
1058  * @param[in] mmax Maximum energy migration (>0).
1059  * @param[in] nmbins Number of migration bins (2 ...).
1060  *
1061  * Sets the nodes for interpolation in migration. None of the nodes will
1062  * have a value of zero unless mmax is zero.
1063  ***************************************************************************/
1064 void GCTACubeEdisp::set_migras(const double& mmax, const int& nmbins)
1065 {
1066  // Clear migration axis
1067  m_migras.clear();
1068 
1069  // Set the migration nodes
1070  double binsize = mmax / double(nmbins);
1071  for (int i = 0; i < nmbins; ++i) {
1072  double migra = binsize * double(i+1); // avoid central singularity
1073  m_migras.append(migra);
1074  }
1075 
1076  // Return
1077  return;
1078 }
1079 
1080 
1081 /***********************************************************************//**
1082  * @brief Compute true energy boundary vector
1083  *
1084  * Computes for all energies of the energy dispersion cube the boundaries in
1085  * true energy that encompass non-zero migration matrix elements. In case
1086  * that no matrix elements are found for a given energy, the interval of true
1087  * energies will be set to [0,0] (i.e. an empty interval).
1088  ***************************************************************************/
1090 {
1091  // Clear energy boundary vector
1092  m_ebounds.clear();
1093 
1094  // Loop over all reconstructed energies
1095  for (int i = 0; i < m_energies.size(); i++) {
1096 
1097  // Initialise results
1098  double migra_min = 0.0;
1099  double migra_max = 0.0;
1100  bool minFound = false;
1101  bool maxFound = false;
1102 
1103  // Compute map of sky ranges to be considered
1104  int mapmin = m_migras.size() * i;
1105  int mapmax = m_migras.size() * (i + 1);
1106 
1107  // Loop over sky map entries from lower migra
1108  for (int imap = mapmin; imap < mapmax; ++imap) {
1109 
1110  // Get maximum energy dispersion term for all sky pixels
1111  double edisp = 0.0;
1112  for (int pixel = 0; pixel < m_cube.npix(); ++pixel) {
1113  double value = m_cube(pixel, imap);
1114  if (value > edisp) {
1115  edisp = value;
1116  }
1117  }
1118 
1119  // Find first non-zero energy dispersion term
1120  if (edisp > 0.0) {
1121  minFound = true;
1122  migra_min = m_migras[imap - mapmin];
1123  break;
1124  }
1125 
1126  } // endfor: loop over maps
1127 
1128  // Loop over sky map entries from high migra
1129  for (int imap = mapmax-1; imap >= mapmin; --imap) {
1130 
1131  // Get maximum energy dispersion term for all sky pixels
1132  double edisp = 0.0;
1133  for (int pixel = 0; pixel < m_cube.npix(); ++pixel) {
1134  double value = m_cube(pixel, imap);
1135  if (value > edisp) {
1136  edisp = value;
1137  }
1138  }
1139 
1140  // Find first non-zero energy dispersion term
1141  if (edisp > 0.0) {
1142  maxFound = true;
1143  migra_max = m_migras[imap - mapmin];
1144  break;
1145  }
1146 
1147  } // endfor: loop over maps
1148 
1149  // Initialise energy boundaries
1150  GEnergy emin;
1151  GEnergy emax;
1152 
1153  // Compute energy boundaries if they were found and if they are
1154  // valid
1155  if (minFound && maxFound && migra_min > 0.0 && migra_max > 0.0 &&
1156  migra_max > migra_min) {
1157  emin = m_energies[i] * migra_min;
1158  emax = m_energies[i] * migra_max;
1159  }
1160 
1161  // ... otherwise we set the interval to a zero interval for safety
1162  else {
1163  emin.clear();
1164  emax.clear();
1165  }
1166 
1167  // Append energy boundaries
1168  m_ebounds.push_back(GEbounds(emin, emax));
1169 
1170  } // endfor: looped over all reconstruced energies
1171 
1172  // Return
1173  return;
1174 }
const GFilename & filename(void) const
Return edisp cube filename.
std::string print(const GChatter &chatter=NORMAL) const
Print energy container information.
Definition: GEnergies.cpp:764
int offset(const int &imigra, const int &iebin) const
Return map offset.
int size(void) const
Return number of nodes in node array.
Definition: GNodeArray.hpp:188
void clear(void)
Clear energy container.
Definition: GEnergies.cpp:225
std::string print(const GChatter &chatter=NORMAL) const
Print energy dispersion cube information.
Sky map class.
Definition: GSkyMap.hpp:89
Abstract FITS image base class.
Definition: GFitsImage.hpp:43
CTA event list class interface definition.
void read(const GFitsTable &table)
Read energies from FITS table.
Definition: GEnergies.cpp:681
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
double norm(const GVector &vector)
Computes vector norm.
Definition: GVector.cpp:821
void write(GFits &file, const std::string &extname=gammalib::extname_energies) const
Write energies into FITS object.
Definition: GEnergies.cpp:725
std::vector< GEbounds > m_ebounds
Energy boundaries.
GFilename m_filename
Filename.
void copy_members(const GCTACubeEdisp &cube)
Copy class members.
void ebounds(const GEbounds &ebounds)
Set energy boundaries.
Definition: GEvents.cpp:136
bool overlaps(const GSkyRegion &region) const
Checks whether a region overlaps with this map.
Definition: GSkyMap.cpp:1944
GEbounds ebounds(void) const
Get energy boundaries.
int m_inx3
Index of upper right node.
void free_members(void)
Delete class members.
CTA instrument response function class definition.
const double & radius(void) const
Returns radius of region of interest in degrees.
Definition: GCTARoi.hpp:124
virtual ~GCTACubeEdisp(void)
Destructor.
const double & wgt_left(void) const
Returns left node weight.
Definition: GNodeArray.hpp:260
int m_inx1
Index of upper left node.
Abstract base class for the CTA energy dispersion.
Definition: GCTAEdisp.hpp:49
void write(GFits &file) const
Write energy dispersion cube into FITS file.
virtual void response(const GResponse &rsp)
Set response function.
void clear(void)
Clear node array.
Definition: GNodeArray.cpp:235
void counts(const GSkyMap &counts)
Set event cube counts from sky map.
double log10TeV(void) const
Return log10 of energy in TeV.
Definition: GEnergy.cpp:457
void fill(const GObservations &obs, GLog *log=NULL)
Fill energy dispersion cube from observation container.
void init_members(void)
Initialise class members.
const GEnergies & energies(void) const
Return energies for energy dispersion cub.
Interface for the CTA region of interest class.
Definition: GCTARoi.hpp:48
void load(const GFilename &filename)
Load energy dispersion cube from FITS file.
Gammalib tools definition.
FITS file class.
Definition: GFits.hpp:63
FITS file class interface definition.
Interface for the circular sky region class.
GSkyMap m_cube
Energy dispersion cube.
void set_value(const double &value) const
Set indices and weighting factors for interpolation.
Definition: GNodeArray.cpp:580
void pointing(const GCTAPointing &pointing)
Set CTA pointing.
GNodeArray m_migras
Migra bins for the Edisp cube.
void read(const GFits &fits)
Read energy dispersion cube from FITS file.
int m_inx4
Index of lower right node.
void read(const GFitsHDU &hdu)
Read skymap from FITS HDU.
Definition: GSkyMap.cpp:2434
const std::vector< int > & shape(void) const
Returns shape of maps.
Definition: GSkyMap.hpp:390
void dir(const GSkyDir &dir)
Set sky direction.
Energy container class.
Definition: GEnergies.hpp:60
#define G_EBOUNDS
void id(const std::string &id)
Set observation identifier.
Abstract FITS image base class definition.
const std::string extname_cta_migras
double migra_max(void) const
Return maximum migra value.
void write(GFits &fits, const std::string &extname="NODES") const
Write nodes into FITS object.
Definition: GNodeArray.cpp:807
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
GCTACubeEdisp(void)
Void constructor.
Information logger interface definition.
Definition: GLog.hpp:62
CTA event bin container class.
virtual double livetime(void) const
Return livetime.
double m_wgt4
Weight of lower right node.
const double & wgt_right(void) const
Returns right node weight.
Definition: GNodeArray.hpp:275
Energy boundaries container class.
Definition: GEbounds.hpp:60
void remove(const int &extno)
Remove HDU from FITS file.
Definition: GFits.cpp:848
GEbounds ebounds(const GEnergy &obsEng) const
Return boundaries in true energy.
const int & nmaps(void) const
Returns number of maps.
Definition: GSkyMap.hpp:378
GSkyDir inx2dir(const int &index) const
Returns sky direction of pixel.
Definition: GSkyMap.cpp:1321
GCTARoi roi(void) const
Get Region of Interest.
void clear_cube(void)
Clear all pixels in the energy dispersion cube.
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
const GCTAInstDir & centre(void) const
Returns region of interest centre.
Definition: GCTARoi.hpp:110
void set(const GCTAObservation &obs)
Set energy dispersion cube for one CTA observation.
int size(void) const
Return number of energies in container.
Definition: GEnergies.hpp:162
virtual std::string instrument(void) const
Return instrument name.
double m_wgt3
Weight of upper right node.
GVector log(const GVector &vector)
Computes natural logarithm of vector elements.
Definition: GVector.cpp:1184
GEnergies m_energies
True energy values of cube.
#define G_CONSTRUCTOR2
Abstract interface for FITS table.
Definition: GFitsTable.hpp:44
CTA energy dispersion for stacked analysis.
GChatter
Definition: GTypemaps.hpp:33
const int & inx_left(void) const
Returns left node index.
Definition: GNodeArray.hpp:231
void clear(void)
Clear instance.
Definition: GSkyMap.cpp:1090
void set_migras(const double &mmax, const int &nmbins)
Set nodes for interpolation in migration.
GCTACubeEdisp & operator=(const GCTACubeEdisp &cube)
Assignment operator.
#define G_CONSTRUCTOR1
double m_wgt1
Weight of upper left node.
int size(void) const
Return number of observations in container.
void compute_ebounds(void) const
Compute true energy boundary vector.
CTA observation class interface definition.
Observation container class.
double operator()(const GEnergy &ereco, const GEnergy &etrue, const GSkyDir &dir) const
Return energy dispersion in units of MeV .
CTA instrument response function class.
const int & inx_right(void) const
Returns right node index.
Definition: GNodeArray.hpp:245
void update(const GEnergy &ereco, const GEnergy &etrue) const
Update energy dispersion cube indices and weights cache.
bool is_empty(void) const
Signals if there are no energies in container.
Definition: GEnergies.hpp:176
const std::string extname_energies
Definition: GEnergies.hpp:44
void set(const GEbounds &ebounds)
Set energies from energy boundaries.
Definition: GEnergies.cpp:416
void name(const std::string &name)
Set observation name.
#define G_FILL_CUBE
GFitsImage * image(const int &extno)
Get pointer to image HDU.
Definition: GFits.cpp:360
Information logger class definition.
void clear(void)
Clear file name.
Definition: GFilename.cpp:188
GCTACubeEdisp * clone(void) const
Clone energy dispersion cube.
GNodeArray m_elogmeans
Mean log10TeV energy for the Edisp cube.
CTA event bin container class interface definition.
void set_eng_axis(void)
Set nodes for interpolation in true energy.
std::string eventtype(void) const
Return event type string.
bool is_valid(void) const
Checks if RoI is valid.
Definition: GCTARoi.hpp:151
double dist(const GSkyDir &dir) const
Compute angular distance between sky directions in radians.
Definition: GSkyDir.hpp:265
double m_wgt2
Weight of lower left node.
void read(const GFitsTable &table)
Read nodes from FITS table.
Definition: GNodeArray.cpp:770
int m_inx2
Index of lower left node.
CTA observation class.
std::string parformat(const std::string &s, const int &indent=0)
Convert string in parameter format.
Definition: GTools.cpp:1022
void clear(void)
Clear energy dispersion cube.
void append(const double &node)
Append one node to array.
Definition: GNodeArray.cpp:351
const int & npix(void) const
Returns number of pixels.
Definition: GSkyMap.hpp:334
void fill_cube(const GCTAObservation &obs, GSkyMap *exposure=NULL, GLog *log=NULL)
Fill energy dispersion cube from observation container.
Observations container class interface definition.
Sky direction class.
Definition: GSkyDir.hpp:62
void close(void)
Close FITS file.
Definition: GFits.cpp:1314
Circular sky region class interface definition.
void save(const GFilename &filename, const bool &clobber=false) const
Save energy dispersion cube into FITS file.
void clear(void)
Clear instance.
Definition: GEnergy.cpp:261
CTA cube analysis energy disperson class definition.
Mathematical function definitions.
Class that handles energies in a unit independent way.
Definition: GEnergy.hpp:48
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.