Stefan Eschbach / gammalib

Commit 804f12cecc6f19874a7d0087a23da2014c911839

Authored by Jürgen Knödlseder
1 parent c5844a9a

Add energy margin in response cube computation

Inline Side-by-side
Showing 4 changed files with 98 additions and 76 deletions
inst/cta/src/GCTACubeBackground.cpp
  View file @804f12c
... ... @@ -57,7 +57,6 @@
57 57 //#define G_DEBUG_CACHE
58 58  
59 59 /* __ Constants __________________________________________________________ */
60   -const GEnergy g_energy_margin(1.0e-12, "TeV");
61 60  
62 61  
63 62 /*==========================================================================
... ... @@ -345,7 +344,8 @@ void GCTACubeBackground::fill(const GObservations& obs, GLog* log)
345 344 for (int i = 0; i < obs.size(); ++i) {
346 345  
347 346 // Get observation and continue only if it is a CTA observation
348   - const GCTAObservation* cta = dynamic_cast<const GCTAObservation*>(obs[i]);
  347 + const GCTAObservation* cta =
  348 + dynamic_cast<const GCTAObservation*>(obs[i]);
349 349  
350 350 // Skip observation if it's not CTA
351 351 if (cta == NULL) {
... ... @@ -726,7 +726,7 @@ void GCTACubeBackground::set_eng_axis(void)
726 726  
727 727 // Compute nodes
728 728 for (int i = 0; i < bins; ++i) {
729   -
  729 +
730 730 // Get logE/TeV
731 731 m_elogmeans.append(m_energies[i].log10TeV());
732 732  
... ...
inst/cta/src/GCTACubeEdisp.cpp
  View file @804f12c
... ... @@ -57,7 +57,6 @@
57 57 /* __ Debug definitions __________________________________________________ */
58 58  
59 59 /* __ Constants __________________________________________________________ */
60   -const GEnergy g_energy_margin(1.0e-12, "TeV");
61 60  
62 61  
63 62 /*==========================================================================
... ... @@ -642,10 +641,10 @@ GEbounds GCTACubeEdisp::ebounds(const GEnergy&amp; obsEng) const
642 641 throw GException::invalid_value(G_EBOUNDS, msg);
643 642 }
644 643  
645   - // If ebounds were not computed, before, compute them now
646   - if (m_ebounds.empty()) {
647   - compute_ebounds();
648   - }
  644 + // If ebounds were not computed, before, compute them now
  645 + if (m_ebounds.empty()) {
  646 + compute_ebounds();
  647 + }
649 648  
650 649 // Return the index of the energy boundaries that are just below the
651 650 // observed energy. As the energy dispersion decreases with increasing
... ... @@ -751,7 +750,7 @@ void GCTACubeEdisp::init_members(void)
751 750 m_wgt3 = 0.0;
752 751 m_wgt4 = 0.0;
753 752 m_ebounds.clear();
754   -
  753 +
755 754 // Return
756 755 return;
757 756 }
... ... @@ -833,6 +832,9 @@ void GCTACubeEdisp::fill_cube(const GCTAObservation&amp; obs,
833 832 GSkyMap* exposure,
834 833 GLog* log)
835 834 {
  835 + // Set energy margin
  836 + static const GEnergy margin(1.0, "MeV");
  837 +
836 838 // Only continue if we have an event list
837 839 if (obs.eventtype() == "EventList") {
838 840  
... ... @@ -906,7 +908,7 @@ void GCTACubeEdisp::fill_cube(const GCTAObservation&amp; obs,
906 908 if (roi.centre().dir().dist_deg(dir) > roi.radius()) {
907 909 continue;
908 910 }
909   -
  911 +
910 912 // Compute theta angle with respect to pointing direction in
911 913 // radians
912 914 double theta = pnt.dist(dir);
... ... @@ -920,7 +922,11 @@ void GCTACubeEdisp::fill_cube(const GCTAObservation&amp; obs,
920 922 // energy boundaries are reconstructed energies, hence
921 923 // this is only an approximation, but probably the only
922 924 // we can really do.
923   - if (!obs_ebounds.contains(m_energies[iebin])) {
  925 + // We allow here for a small margin in case of rounding
  926 + // errors in the energy boundaries.
  927 + if (!(obs_ebounds.contains(m_energies[iebin]) ||
  928 + obs_ebounds.contains(m_energies[iebin]-margin) ||
  929 + obs_ebounds.contains(m_energies[iebin]+margin))) {
924 930 continue;
925 931 }
926 932  
... ... @@ -961,7 +967,7 @@ void GCTACubeEdisp::fill_cube(const GCTAObservation&amp; obs,
961 967 logEsrc) * weight;
962 968  
963 969 } // endfor: looped over migration bins
964   -
  970 +
965 971 } // endfor: looped over energy bins
966 972  
967 973 } // endfor: looped over all pixels
... ... @@ -989,9 +995,9 @@ void GCTACubeEdisp::update(const double&amp; migra, const double&amp; logEsrc) const
989 995 {
990 996 // Set node array for energy interpolation
991 997 m_elogmeans.set_value(logEsrc);
992   -
  998 +
993 999 // Set node array for delta interpolation
994   - m_migras.set_value(migra);
  1000 + m_migras.set_value(migra);
995 1001  
996 1002 // Set indices for bi-linear interpolation
997 1003 m_inx1 = offset(m_migras.inx_left(), m_elogmeans.inx_left());
... ... @@ -1022,7 +1028,7 @@ void GCTACubeEdisp::set_eng_axis(void)
1022 1028 for (int i = 0; i < m_energies.size(); ++i) {
1023 1029 m_elogmeans.append(m_energies[i].log10TeV());
1024 1030 }
1025   -
  1031 +
1026 1032 // Return
1027 1033 return;
1028 1034 }
... ... @@ -1039,29 +1045,29 @@ void GCTACubeEdisp::set_eng_axis(void)
1039 1045 ***************************************************************************/
1040 1046 void GCTACubeEdisp::compute_ebounds() const
1041 1047 {
1042   - // Set epsilon
1043   - const double eps = 1.0e-12;
  1048 + // Set epsilon
  1049 + const double eps = 1.0e-12;
1044 1050  
1045 1051 // Clear energy boundary vector
1046 1052 m_ebounds.clear();
1047 1053  
1048 1054 // Loop over all energies
1049   - for (int i = 0; i < m_energies.size(); i++) {
  1055 + for (int i = 0; i < m_energies.size(); i++) {
1050 1056  
1051   - // Initialise results
1052   - double migra_min = 0.0;
1053   - double migra_max = 0.0;
1054   - bool minFound = false;
1055   - bool maxFound = false;
  1057 + // Initialise results
  1058 + double migra_min = 0.0;
  1059 + double migra_max = 0.0;
  1060 + bool minFound = false;
  1061 + bool maxFound = false;
1056 1062  
1057   - // Compute map of sky ranges to be considered
1058   - int mapmin = m_migras.size() * i;
1059   - int mapmax = m_migras.size() * (i + 1);
  1063 + // Compute map of sky ranges to be considered
  1064 + int mapmin = m_migras.size() * i;
  1065 + int mapmax = m_migras.size() * (i + 1);
1060 1066  
1061   - // Loop over sky map entries from lower migra
1062   - for (int imap = mapmin; imap < mapmax; ++imap) {
  1067 + // Loop over sky map entries from lower migra
  1068 + for (int imap = mapmin; imap < mapmax; ++imap) {
1063 1069  
1064   - // Get maximum energy dispersion term for all sky pixels
  1070 + // Get maximum energy dispersion term for all sky pixels
1065 1071 double edisp = 0.0;
1066 1072 for (int pixel = 0; pixel < m_cube.npix(); ++pixel) {
1067 1073 double value = m_cube(pixel, imap);
... ... @@ -1070,19 +1076,19 @@ void GCTACubeEdisp::compute_ebounds() const
1070 1076 }
1071 1077 }
1072 1078  
1073   - // Find first non-negligible energy dispersion term
1074   - if (edisp >= eps) {
1075   - minFound = true;
1076   - migra_min = m_migras[imap - mapmin];
1077   - break;
1078   - }
  1079 + // Find first non-negligible energy dispersion term
  1080 + if (edisp >= eps) {
  1081 + minFound = true;
  1082 + migra_min = m_migras[imap - mapmin];
  1083 + break;
  1084 + }
1079 1085  
1080   - } // endfor: loop over maps
  1086 + } // endfor: loop over maps
1081 1087  
1082   - // Loop over sky map entries from high migra
1083   - for (int imap = mapmax-1; imap >= mapmin; --imap) {
  1088 + // Loop over sky map entries from high migra
  1089 + for (int imap = mapmax-1; imap >= mapmin; --imap) {
1084 1090  
1085   - // Get maximum energy dispersion term for all sky pixels
  1091 + // Get maximum energy dispersion term for all sky pixels
1086 1092 double edisp = 0.0;
1087 1093 for (int pixel = 0; pixel < m_cube.npix(); ++pixel) {
1088 1094 double value = m_cube(pixel, imap);
... ... @@ -1091,37 +1097,37 @@ void GCTACubeEdisp::compute_ebounds() const
1091 1097 }
1092 1098 }
1093 1099  
1094   - // Find first non-negligible energy dispersion term
1095   - if (edisp >= eps) {
1096   - maxFound = true;
1097   - migra_max = m_migras[imap - mapmin];
1098   - break;
1099   - }
  1100 + // Find first non-negligible energy dispersion term
  1101 + if (edisp >= eps) {
  1102 + maxFound = true;
  1103 + migra_max = m_migras[imap - mapmin];
  1104 + break;
  1105 + }
1100 1106  
1101   - } // endfor: loop over maps
  1107 + } // endfor: loop over maps
1102 1108  
1103   - // Initialise energy boundaries
1104   - GEnergy emin;
1105   - GEnergy emax;
  1109 + // Initialise energy boundaries
  1110 + GEnergy emin;
  1111 + GEnergy emax;
1106 1112  
1107   - // Compute energy boundaries if they were found and if they are
  1113 + // Compute energy boundaries if they were found and if they are
1108 1114 // valid
1109   - if (minFound && maxFound && migra_min > 0.0 && migra_max > 0.0 &&
  1115 + if (minFound && maxFound && migra_min > 0.0 && migra_max > 0.0 &&
1110 1116 migra_max > migra_min) {
1111   - emin = m_energies[i] * migra_min;
1112   - emax = m_energies[i] * migra_max;
1113   - }
1114   -
  1117 + emin = m_energies[i] * migra_min;
  1118 + emax = m_energies[i] * migra_max;
  1119 + }
  1120 +
1115 1121 // ... otherwise we set the interval to a zero interval for safety
1116   - else {
1117   - emin.log10TeV(0.0);
1118   - emax.log10TeV(0.0);
1119   - }
  1122 + else {
  1123 + emin.log10TeV(0.0);
  1124 + emax.log10TeV(0.0);
  1125 + }
1120 1126  
1121   - // Append energy boundaries
1122   - m_ebounds.push_back(GEbounds(emin, emax));
  1127 + // Append energy boundaries
  1128 + m_ebounds.push_back(GEbounds(emin, emax));
1123 1129  
1124   - } // endfor: looped over all energies
  1130 + } // endfor: looped over all energies
1125 1131  
1126 1132 // Return
1127 1133 return;
... ...
inst/cta/src/GCTACubeExposure.cpp
  View file @804f12c
... ... @@ -238,7 +238,7 @@ GCTACubeExposure&amp; GCTACubeExposure::operator=(const GCTACubeExposure&amp; cube)
238 238 * @return Exposure (in units of cm2 s)
239 239 ***************************************************************************/
240 240 double GCTACubeExposure::operator()(const GSkyDir& dir, const GEnergy& energy) const
241   -{
  241 +{
242 242 // Set indices and weighting factors for interpolation
243 243 update(energy.log10TeV());
244 244  
... ... @@ -364,7 +364,7 @@ void GCTACubeExposure::fill(const GObservations&amp; obs, GLog* log)
364 364  
365 365 // Fill exposure cube
366 366 fill_cube(*cta, log);
367   -
  367 +
368 368 } // endfor: looped over observations
369 369  
370 370 // Return
... ... @@ -566,7 +566,7 @@ void GCTACubeExposure::init_members(void)
566 566 m_inx_right = 0;
567 567 m_wgt_left = 0.0;
568 568 m_wgt_right = 0.0;
569   -
  569 +
570 570 // Return
571 571 return;
572 572 }
... ... @@ -621,6 +621,9 @@ void GCTACubeExposure::free_members(void)
621 621 ***************************************************************************/
622 622 void GCTACubeExposure::fill_cube(const GCTAObservation& obs, GLog* log)
623 623 {
  624 + // Set energy margin
  625 + static const GEnergy margin(1.0, "MeV");
  626 +
624 627 // Only continue if we have an event list
625 628 if (obs.eventtype() == "EventList") {
626 629  
... ... @@ -683,7 +686,7 @@ void GCTACubeExposure::fill_cube(const GCTAObservation&amp; obs, GLog* log)
683 686 if (roi.centre().dir().dist_deg(dir) > roi.radius()) {
684 687 continue;
685 688 }
686   -
  689 +
687 690 // Compute theta angle with respect to pointing direction in
688 691 // radians
689 692 double theta = pnt.dist(dir);
... ... @@ -696,7 +699,13 @@ void GCTACubeExposure::fill_cube(const GCTAObservation&amp; obs, GLog* log)
696 699 // energies while the observation energy boundaries are
697 700 // reconstructed energies, hence this is only an
698 701 // approximation, but probably the only we can really do.
699   - if (!obs_ebounds.contains(m_energies[iebin])) {
  702 + // We allow here for a small margin in case of rounding
  703 + // errors in the energy boundaries.
  704 +//std::cout << obs_ebounds.emin() << " <= " << m_energies[iebin] << " <= " << obs_ebounds.emax() << std::endl;
  705 + if (!(obs_ebounds.contains(m_energies[iebin]) ||
  706 + obs_ebounds.contains(m_energies[iebin]-margin) ||
  707 + obs_ebounds.contains(m_energies[iebin]+margin))) {
  708 +//std::cout << "Skip" << std::endl;
700 709 continue;
701 710 }
702 711  
... ... @@ -766,7 +775,7 @@ void GCTACubeExposure::set_eng_axis(void)
766 775  
767 776 // Compute nodes
768 777 for (int i = 0; i < bins; ++i) {
769   -
  778 +
770 779 // Get logE/TeV
771 780 m_elogmeans.append(m_energies[i].log10TeV());
772 781  
... ...
inst/cta/src/GCTACubePsf.cpp
  View file @804f12c
... ... @@ -267,7 +267,7 @@ GCTACubePsf::~GCTACubePsf(void)
267 267 * @param[in] cube Mean PSF cube.
268 268 * @return Mean PSF cube.
269 269 ***************************************************************************/
270   -GCTACubePsf& GCTACubePsf::operator= (const GCTACubePsf& cube)
  270 +GCTACubePsf& GCTACubePsf::operator=(const GCTACubePsf& cube)
271 271 {
272 272 // Execute only if object is not identical
273 273 if (this != &cube) {
... ... @@ -669,7 +669,7 @@ void GCTACubePsf::init_members(void)
669 669 m_wgt2 = 0.0;
670 670 m_wgt3 = 0.0;
671 671 m_wgt4 = 0.0;
672   -
  672 +
673 673 // Return
674 674 return;
675 675 }
... ... @@ -752,6 +752,9 @@ void GCTACubePsf::fill_cube(const GCTAObservation&amp; obs,
752 752 GSkyMap* exposure,
753 753 GLog* log)
754 754 {
  755 + // Set energy margin
  756 + static const GEnergy margin(1.0, "MeV");
  757 +
755 758 // Only continue if we have an event list
756 759 if (obs.eventtype() == "EventList") {
757 760  
... ... @@ -809,7 +812,7 @@ void GCTACubePsf::fill_cube(const GCTAObservation&amp; obs,
809 812  
810 813 // Get pixel sky direction
811 814 GSkyDir dir = m_cube.inx2dir(pixel);
812   -
  815 +
813 816 // Skip pixel if it is outside the RoI
814 817 if (roi.centre().dir().dist_deg(dir) > roi.radius()) {
815 818 continue;
... ... @@ -827,7 +830,11 @@ void GCTACubePsf::fill_cube(const GCTAObservation&amp; obs,
827 830 // energies while the observation energy boundaries are
828 831 // reconstructed energies, hence this is only an
829 832 // approximation, but probably the only we can really do.
830   - if (!obs_ebounds.contains(m_energies[iebin])) {
  833 + // We allow here for a small margin in case of rounding
  834 + // errors in the energy boundaries.
  835 + if (!(obs_ebounds.contains(m_energies[iebin]) ||
  836 + obs_ebounds.contains(m_energies[iebin]-margin) ||
  837 + obs_ebounds.contains(m_energies[iebin]+margin))) {
831 838 continue;
832 839 }
833 840  
... ... @@ -876,7 +883,7 @@ void GCTACubePsf::fill_cube(const GCTAObservation&amp; obs,
876 883 *
877 884 * @param[in] delta Angular separation between true and measured photon
878 885 * directions (radians).
879   - * @param[in] logE Log10 true photon energy (TeV).
  886 + * @param[in] logE Log10 true photon energy (TeV).
880 887 *
881 888 * This method updates the PSF parameter cache.
882 889 ***************************************************************************/
... ... @@ -892,7 +899,7 @@ void GCTACubePsf::update(const double&amp; delta, const double&amp; logE) const
892 899  
893 900 // Set node array for energy interpolation
894 901 m_elogmeans.set_value(logE);
895   -
  902 +
896 903 // Set indices for bi-linear interpolation
897 904 m_inx1 = offset(m_deltas_cache.inx_left(), m_elogmeans.inx_left());
898 905 m_inx2 = offset(m_deltas_cache.inx_left(), m_elogmeans.inx_right());
... ... @@ -986,7 +993,7 @@ void GCTACubePsf::set_eng_axis(void)
986 993  
987 994 // Compute nodes
988 995 for (int i = 0; i < bins; ++i) {
989   -
  996 +
990 997 // Get logE/TeV
991 998 m_elogmeans.append(m_energies[i].log10TeV());
992 999  
... ... @@ -1016,7 +1023,7 @@ void GCTACubePsf::set_to_smooth(void)
1016 1023 m_cube(pixel, idst) = m_cube(pixel, isrc);
1017 1024 }
1018 1025 }
1019   -
  1026 +
1020 1027 // Get index of last delta bin
1021 1028 int idelta = m_deltas.size()-1;
1022 1029  
... ...