Justine Devin / gammalib

Commit 528cae1913fcec4dd6985df6aa6560b24156aa70

Authored by Jürgen Knödlseder
1 parent e1b87df5

Correct GCTAEdisp2D::compute_ebounds_obs() and GCTAEdisp2D::compute_ebouds_src()… 

… method so that they return the correctly return the energy intervals covered by valid matrix elements (#1662).

The former code was not limiting the maximum to the elements covered by the matrix, leading eventually to integration warnings.
Inline Side-by-side
Showing 3 changed files with 96 additions and 50 deletions
ChangeLog
  View file @528cae1
1   -2016-01-18 Juergen Knoedlseder <jurgen.knodlseder@irap.omp.eu>
  1 +2016-02-22 Juergen Knoedlseder <jurgen.knodlseder@irap.omp.eu>
2 2  
3 3 * Version 1.1.0 released
4 4 =========================
5 5  
  6 + Correct GCTAEdisp2D::compute_ebounds methods (#1662)
6 7 Use GFilename in GFits to store filename (#1649)
7 8 Load CTA GTI extension from the HDU name provided by the DSS keyword (#1598)
8 9  
... ...
... ... @@ -2,7 +2,7 @@ New Features and Important Changes in GammaLib 1.1.0
2 2  
3 3 Author(s) : Juergen Knoedlseder
4 4  
5   -20 January 2016
  5 +11 February 2016
6 6  
7 7  
8 8 1. Introduction
... ... @@ -123,6 +123,11 @@ optionally specified extension in the format [EVENTS;GTI] into an events and a
123 123 Good Time Intervals extension name. These changes satisfy the feature request
124 124 #1598.
125 125  
  126 +The GCTAEdisp2D::compute_ebounds_obs() and GCTAEdisp2::compute_ebounds_src()
  127 +methods have been corrected so that they return the energy intervals covered by
  128 +valid matrix elements. With the code before there were problems with not setting
  129 +correctly the maximum energies of the intervals (#1662).
  130 +
126 131  
127 132 21. COMPTEL interface
128 133 ---------------------
... ...
inst/cta/src/GCTAEdisp2D.cpp
  View file @528cae1
... ... @@ -725,6 +725,12 @@ void GCTAEdisp2D::free_members(void)
725 725 * @param[in] phi Azimuth angle (rad).
726 726 * @param[in] zenith Zenith angle (rad).
727 727 * @param[in] azimuth Azimuth angle (rad).
  728 + *
  729 + * Computes for all true energies the energy boundaries of the observed
  730 + * energies covered by valid migration matrix elements. Only matrix elements
  731 + * with values >= 1.0e-12 are considered as valid elements. In case that no
  732 + * matrix elements are found of a given true energy, the interval of observed
  733 + * energies will be set to [1 TeV, 1 TeV] (i.e. an empty interval).
728 734 ***************************************************************************/
729 735 void GCTAEdisp2D::compute_ebounds_obs(const double& theta,
730 736 const double& phi,
... ... @@ -741,49 +747,62 @@ void GCTAEdisp2D::compute_ebounds_obs(const double&amp; theta,
741 747 for (int isrc = 0; isrc < m_edisp.axis(0); ++isrc) {
742 748  
743 749 // Set Esrc
744   - double Esrc = std::sqrt(m_edisp.axis_hi(0, isrc) *
745   - m_edisp.axis_lo(0, isrc));
  750 + double Esrc = std::sqrt(m_edisp.axis_hi(0,isrc) *
  751 + m_edisp.axis_lo(0,isrc));
746 752 double logEsrc = std::log10(Esrc);
747 753  
748 754 // Initialise results
749   - double logEobsMin = -30.0;
750   - double logEobsMax = 10.0;
  755 + double logEobsMin = 0.0;
  756 + double logEobsMax = 0.0;
751 757 bool minFound = false;
752 758 bool maxFound = false;
753 759  
754   - // Find boundaries, loop over EobsOverEsrc
755   - for (int i = 0; i < m_edisp.axis(1); ++i) {
  760 + // Determine number of MIGRA bins
  761 + int n_migra = m_edisp.axis(1);
  762 +
  763 + // Find minimum boundary
  764 + for (int i = 0; i < n_migra; ++i) {
  765 +
  766 + // Compute EobsOverEsrc value
  767 + double EobsOverEsrc = 0.5 * (m_edisp.axis_hi(1,i) +
  768 + m_edisp.axis_lo(1,i));
756 769  
757   - // Compute value EobsOverEsrc
758   - double EobsOverEsrc = 0.5 * (m_edisp.axis_hi(1, i) +
759   - m_edisp.axis_lo(1, i));
  770 + // Get matrix term
  771 + double edisp = m_edisp(0, logEsrc, EobsOverEsrc, theta);
760 772  
761 773 // Find first non-negligible matrix term
762   - if (!minFound && m_edisp(0, logEsrc, EobsOverEsrc, theta) >= eps) {
  774 + if (edisp >= eps) {
763 775 minFound = true;
764 776 logEobsMin = std::log10(EobsOverEsrc * Esrc);
  777 + break;
765 778 }
766 779  
767   - // Find last non-negligible matrix term
768   - else if (minFound && !maxFound &&
769   - m_edisp(0, logEsrc, EobsOverEsrc, theta) < eps) {
  780 + } // endfor: find minimum boundary
  781 +
  782 + // Find maximum boundary
  783 + for (int i = n_migra-1; i >= 0; i--) {
  784 +
  785 + // Compute EobsOverEsrc value
  786 + double EobsOverEsrc = 0.5 * (m_edisp.axis_hi(1,i) +
  787 + m_edisp.axis_lo(1,i));
  788 +
  789 + // Get matrix term
  790 + double edisp = m_edisp(0, logEsrc, EobsOverEsrc, theta);
  791 +
  792 + // Find first non-negligible matrix term
  793 + if (edisp >= eps) {
770 794 maxFound = true;
771 795 logEobsMax = std::log10(EobsOverEsrc * Esrc);
  796 + break;
772 797 }
773 798  
774   - // Continue searching
775   - else if (minFound && maxFound &&
776   - m_edisp(0, logEsrc, EobsOverEsrc, theta) >= eps) {
777   - maxFound = false;
778   - }
779   -
780   - } // endfor: looped over EobsOverEsrc
781   -
  799 + } // endfor: find maximum boundary
782 800  
783   - // If energy dispersion has never become negligible until end of loop,
784   - // reset logEobsMax
785   - if (!maxFound) {
786   - logEobsMax = 10.0;
  801 + // If we did not find boundaries then set the interval to
  802 + // a zero interval for safety
  803 + if (!minFound || !maxFound) {
  804 + logEobsMin = 0.0;
  805 + logEobsMax = 0.0;
787 806 }
788 807  
789 808 // Set energy boundaries
... ... @@ -809,6 +828,12 @@ void GCTAEdisp2D::compute_ebounds_obs(const double&amp; theta,
809 828 * @param[in] phi Azimuth angle (rad).
810 829 * @param[in] zenith Zenith angle (rad).
811 830 * @param[in] azimuth Azimuth angle (rad).
  831 + *
  832 + * Computes for all observed energies the energy boundaries of the true
  833 + * energies covered by valid migration matrix elements. Only matrix elements
  834 + * with values >= 1.0e-12 are considered as valid elements. In case that no
  835 + * matrix elements are found of a given observed energy, the interval of true
  836 + * energies will be set to [1 TeV, 1 TeV] (i.e. an empty interval).
812 837 ***************************************************************************/
813 838 void GCTAEdisp2D::compute_ebounds_src(const double& theta,
814 839 const double& phi,
... ... @@ -825,49 +850,64 @@ void GCTAEdisp2D::compute_ebounds_src(const double&amp; theta,
825 850 for (int iobs = 0; iobs < m_edisp.axis(0); ++iobs) {
826 851  
827 852 // Set Eobs
828   - double Eobs = std::sqrt(m_edisp.axis_hi(0, iobs) *
829   - m_edisp.axis_lo(0, iobs));
  853 + double Eobs = std::sqrt(m_edisp.axis_hi(0,iobs) *
  854 + m_edisp.axis_lo(0,iobs));
830 855 double logEobs = std::log10(Eobs);
831 856  
832 857 // Initialise results
833   - double logEsrcMin = -10.0;
834   - double logEsrcMax = 30.0;
  858 + double logEsrcMin = 0.0;
  859 + double logEsrcMax = 0.0;
835 860 bool minFound = false;
836 861 bool maxFound = false;
837 862  
838   - // Find boundaries, loop over Esrc
839   - for (int isrc = 0; isrc < m_edisp.axis(0); ++isrc) {
  863 + // Determine number of true energy bins
  864 + int n_bins = m_edisp.axis(0);
  865 +
  866 + // Find minimum boundary
  867 + for (int isrc = 0; isrc < n_bins; ++isrc) {
840 868  
841 869 // Set Esrc
842   - double Esrc = std::sqrt(m_edisp.axis_hi(0, isrc) *
843   - m_edisp.axis_lo(0, isrc));
  870 + double Esrc = std::sqrt(m_edisp.axis_hi(0,isrc) *
  871 + m_edisp.axis_lo(0,isrc));
844 872 double logEsrc = std::log10(Esrc);
845 873  
  874 + // Get matrix term
  875 + double edisp = operator()(logEobs, logEsrc, theta);
  876 +
846 877 // Find first non-negligible matrix term
847   - if (!minFound && operator()(logEobs, logEsrc, theta) >= eps) {
  878 + if (edisp >= eps) {
848 879 minFound = true;
849 880 logEsrcMin = logEsrc;
  881 + break;
850 882 }
851 883  
852   - // Find last non-negligible matrix term
853   - else if (minFound && !maxFound &&
854   - operator()(logEobs, logEsrc, theta) < eps) {
  884 + } // endfor: find minimum boundary
  885 +
  886 + // Find maximum boundary
  887 + for (int isrc = n_bins-1; isrc >= 0; isrc--) {
  888 +
  889 + // Set Esrc
  890 + double Esrc = std::sqrt(m_edisp.axis_hi(0,isrc) *
  891 + m_edisp.axis_lo(0,isrc));
  892 + double logEsrc = std::log10(Esrc);
  893 +
  894 + // Get matrix term
  895 + double edisp = operator()(logEobs, logEsrc, theta);
  896 +
  897 + // Find first non-negligible matrix term
  898 + if (edisp >= eps) {
855 899 maxFound = true;
856 900 logEsrcMax = logEsrc;
  901 + break;
857 902 }
858 903  
859   - // Continue
860   - else if (minFound && maxFound &&
861   - operator()(logEobs, logEsrc, theta) >= eps) {
862   - maxFound = false;
863   - }
864   -
865   - } // endfor: looped over true energy
  904 + } // endfor: find maximum boundary
866 905  
867   - // If energy dispersion has never become negligible until end of loop,
868   - // reset logEobsMax
869   - if (!maxFound) {
870   - logEsrcMax = 30.0;
  906 + // If we did not find boundaries then set the interval to
  907 + // a zero interval for safety
  908 + if (!minFound || !maxFound) {
  909 + logEsrcMin = 0.0;
  910 + logEsrcMax = 0.0;
871 911 }
872 912  
873 913 // Set energy boundaries
... ...