Simon Bonnefoy / gammalib

Commit cee203afc123c88a70407ce6785bf815c3a6ddd3

Authored by Florent Forest
1 parent aa568ec0

Implemented rejection in GCTAEdisp2D and GCTAEdispRmf mc methods 

Implemented and modified ebounds_src and ebounds_obs methods in GCTAEdisp2D and GCTAEdispRmf
Corrected edisp_matrix.fits file
Inline Side-by-side
Showing 11 changed files with 827 additions and 271 deletions
include/GRmf.hpp
  View file @cee203a
... ... @@ -81,6 +81,8 @@ public:
81 81 const std::string& unit = "keV") const;
82 82 const std::string& filename(void) const;
83 83 std::string print(const GChatter& chatter = NORMAL) const;
  84 + int itruemax(void) const;
  85 + int imeasmax(void) const;
84 86  
85 87 protected:
86 88 // Protected methods
... ... @@ -93,6 +95,8 @@ protected:
93 95 GEbounds m_ebds_true; //!< True energy boundaries
94 96 GEbounds m_ebds_measured; //!< Measured energy boundaries
95 97 GMatrixSparse m_matrix; //!< Sparse redistribution matrix
  98 + mutable int m_imeasmax;
  99 + mutable int m_itruemax;
96 100  
97 101 };
98 102  
... ... @@ -239,4 +243,31 @@ const std::string&amp; GRmf::filename(void) const
239 243 return (m_filename);
240 244 }
241 245  
  246 +
  247 +/***********************************************************************//**
  248 + * @brief Return true energy index of maximum value of the redistribution matrix
  249 + *
  250 + * @return True energy index of maximum value of the redistribution matrix.
  251 + *
  252 + * Returns the true energy index of maximum value of the redistribution matrix.
  253 + ***************************************************************************/
  254 +inline
  255 +int GRmf::itruemax(void) const
  256 +{
  257 + return m_itruemax;
  258 +}
  259 +
  260 +/***********************************************************************//**
  261 + * @brief Return measured energy index of maximum value of the redistribution matrix
  262 + *
  263 + * @return Measured energy index of maximum value of the redistribution matrix.
  264 + *
  265 + * Returns the measured energy index of maximum value of the redistribution matrix.
  266 + ***************************************************************************/
  267 +inline
  268 +int GRmf::imeasmax(void) const
  269 +{
  270 + return m_imeasmax;
  271 +}
  272 +
242 273 #endif /* GRMF_HPP */
... ...
inst/cta/include/GCTAEdisp2D.hpp
  View file @cee203a
... ... @@ -107,6 +107,15 @@ private:
107 107 const double& phi = 0.0,
108 108 const double& zenith = 0.0,
109 109 const double& azimuth = 0.0) const;
  110 + void compute_ebounds_obs(const double& theta = 0.0,
  111 + const double& phi = 0.0,
  112 + const double& zenith = 0.0,
  113 + const double& azimuth = 0.0) const;
  114 + void compute_ebounds_src(const double& theta = 0.0,
  115 + const double& phi = 0.0,
  116 + const double& zenith = 0.0,
  117 + const double& azimuth = 0.0) const;
  118 + void set_fmax(void) const;
110 119  
111 120 // Members
112 121 std::string m_filename; //!< Name of Edisp response file
... ... @@ -114,7 +123,19 @@ private:
114 123  
115 124 mutable bool m_cdf_computed;
116 125 mutable double m_theta;
  126 + mutable bool m_ebounds_obs_computed;
  127 + mutable bool m_ebounds_src_computed;
  128 + mutable double m_logEsrc;
  129 + mutable double m_logEobs;
  130 + mutable int m_index_obs;
  131 + mutable int m_index_src;
  132 + mutable double m_fmax;
  133 +
  134 + mutable std::vector<GEbounds> m_ebounds_obs;
  135 + mutable std::vector<GEbounds> m_ebounds_src;
117 136 mutable std::vector<std::vector<std::pair<double, double> > > m_cumul;
  137 + mutable std::vector<double> m_eobs_axis;
  138 + mutable std::vector<std::vector<double> > m_temp;
118 139 };
119 140  
120 141  
... ...
inst/cta/include/GCTAEdispRmf.hpp
  View file @cee203a
... ... @@ -104,6 +104,14 @@ private:
104 104 const double& phi = 0.0,
105 105 const double& zenith = 0.0,
106 106 const double& azimuth = 0.0) const;
  107 + void compute_ebounds_obs(const double& theta = 0.0,
  108 + const double& phi = 0.0,
  109 + const double& zenith = 0.0,
  110 + const double& azimuth = 0.0) const;
  111 + void compute_ebounds_src(const double& theta = 0.0,
  112 + const double& phi = 0.0,
  113 + const double& zenith = 0.0,
  114 + const double& azimuth = 0.0) const;
107 115  
108 116 // Members
109 117 std::string m_filename; //!< Name of response file
... ... @@ -129,8 +137,17 @@ private:
129 137 mutable std::vector<GVector> m_mc_measured_cdf;
130 138 mutable bool m_cdf_computed;
131 139 mutable double m_theta;
  140 + mutable double m_logEsrc;
  141 + mutable double m_logEobs;
132 142 mutable std::vector<std::vector<std::pair<double, double> > > m_cumul;
133 143  
  144 + mutable bool m_ebounds_obs_computed;
  145 + mutable std::vector<GEbounds> m_ebounds_obs;
  146 + mutable int m_index_obs;
  147 + mutable bool m_ebounds_src_computed;
  148 + mutable std::vector<GEbounds> m_ebounds_src;
  149 + mutable int m_index_src;
  150 +
134 151 };
135 152  
136 153  
... ...
inst/cta/src/GCTAEdisp2D.cpp
  View file @cee203a
... ... @@ -37,6 +37,9 @@
37 37 #include "GFitsBinTable.hpp"
38 38 #include "GCTAEdisp2D.hpp"
39 39  
  40 +//#define G_CDF_METHOD
  41 +#define G_REJECTION
  42 +
40 43 /* __ Method name definitions ____________________________________________ */
41 44  
42 45 /* __ Macros _____________________________________________________________ */
... ... @@ -183,7 +186,7 @@ double GCTAEdisp2D::operator()(const double&amp; logEobs,
183 186 double EobsOverEsrc = std::exp((logEobs-logEsrc) * gammalib::ln10);
184 187  
185 188 // Compute edisp
186   - edisp = m_edisp(0, logEsrc, EobsOverEsrc, theta) * gammalib::ln10;
  189 + edisp = m_edisp(0, logEsrc, EobsOverEsrc, theta);
187 190  
188 191 // Return
189 192 return edisp;
... ... @@ -321,12 +324,29 @@ void GCTAEdisp2D::read(const GFits&amp; fits)
321 324 // Set true energy axis to logarithmic scale
322 325 m_edisp.axis_log10(0);
323 326  
324   - // Set emigra axis to log scale
325   - m_edisp.axis_log10(1);
326   -
327 327 // Set offset angle axis to radians
328 328 m_edisp.axis_radians(2);
329 329  
  330 +
  331 +// for (int i_theta = 0; i_theta < theta_size; ++i_theta) {
  332 +
  333 + m_temp.clear();
  334 +
  335 + int step = migra_size / etrue_size;
  336 +
  337 + for (int itrue = 0; itrue < etrue_size; ++itrue) {
  338 + double logEsrc = std::log10(0.5 * (m_edisp.axis_hi(0, itrue) +
  339 + m_edisp.axis_lo(0, itrue)));
  340 + m_temp.push_back(std::vector<double>());
  341 + for (int iobs = 0; iobs < etrue_size; ++iobs) {
  342 + double EobsOverEsrc = 0.5 * (m_edisp.axis_hi(1, iobs*step) +
  343 + m_edisp.axis_lo(1, iobs*step));
  344 + m_temp[itrue].push_back(operator()(logEsrc, EobsOverEsrc));
  345 + }
  346 + }
  347 +
  348 + // }
  349 +
330 350 // Return
331 351 return;
332 352 }
... ... @@ -375,6 +395,56 @@ void GCTAEdisp2D::save(const std::string&amp; filename, const bool&amp; clobber) const
375 395 return;
376 396 }
377 397  
  398 +/***********************************************************************//**
  399 + * @brief Find maximum energy dispersion
  400 + *
  401 +***************************************************************************/
  402 +void GCTAEdisp2D::set_fmax(void) const
  403 +{
  404 +
  405 + double max = 0.0;
  406 +
  407 + // Loop over Esrc
  408 + for (int i = 0; i < m_edisp.axis(0); ++i) {
  409 +
  410 + // Set Esrc
  411 + double Esrc = 0.5 * (m_edisp.axis_hi(0, i) +
  412 + m_edisp.axis_lo(0, i));
  413 + double logEsrc = std::log10(Esrc);
  414 +
  415 + // Loop over EobsOverEsrc
  416 + for (int j = 0; j < m_edisp.axis(1); ++j) {
  417 +
  418 + // Set EobsOverEsrc
  419 + double EobsOverEsrc = 0.5 * (m_edisp.axis_hi(1, j) +
  420 + m_edisp.axis_lo(1, j));
  421 +/*
  422 + // Loop over theta
  423 + for (int k = 0; k < m_edisp.axis(2); ++k) {
  424 +
  425 + // Set theta
  426 + double theta = 0.5 * (m_edisp.axis_hi(2, k) +
  427 + m_edisp.axis_lo(2, k));
  428 +
  429 + double test = operator()(logEsrc, EobsOverEsrc, theta);
  430 +*/
  431 + double test = operator()(logEsrc, EobsOverEsrc);
  432 + if (max < test) {
  433 + max = test;
  434 + }
  435 +/*
  436 + } // endfor: looped over theta
  437 +*/
  438 + } // endfor: looped over EobsOverEsrc
  439 +
  440 + } // endfor: looped over Esrc
  441 +
  442 + // Set fmax
  443 + m_fmax = max;
  444 +
  445 + // Return
  446 + return;
  447 +}
378 448  
379 449 /***********************************************************************//**
380 450 * @brief Simulate energy dispersion
... ... @@ -388,6 +458,49 @@ void GCTAEdisp2D::save(const std::string&amp; filename, const bool&amp; clobber) const
388 458 *
389 459 * Draws observed energy value.
390 460 ***************************************************************************/
  461 +#if defined(G_REJECTION)
  462 +GEnergy GCTAEdisp2D::mc(GRan& ran,
  463 + const double& logEsrc,
  464 + const double& theta,
  465 + const double& phi,
  466 + const double& zenith,
  467 + const double& azimuth) const
  468 +{
  469 +
  470 + // Set true energy
  471 + GEnergy energy;
  472 +
  473 + // Get boundaries for rejection method
  474 + double emin = ebounds_obs(logEsrc, theta, phi, zenith, azimuth).emin().log10TeV();
  475 +
  476 + double emax = ebounds_obs(logEsrc, theta, phi, zenith, azimuth).emax().log10TeV();
  477 +
  478 + // Set fmax if it has not already been set
  479 + if (m_fmax == -1.0) {
  480 + set_fmax();
  481 + }
  482 +
  483 + double fmax = m_fmax;
  484 +
  485 + // Find energy by rejection method
  486 + double ewidth = emax - emin;
  487 + double e = emin;
  488 + double f = 0.0;
  489 + double ftest = 1.0;
  490 + while (ftest > f) {
  491 + e = emin + ewidth * ran.uniform();
  492 + f = operator()(e, logEsrc, theta, phi, zenith, azimuth);
  493 + ftest = ran.uniform() * fmax;
  494 + }
  495 +
  496 + energy.log10TeV(e);
  497 +
  498 +
  499 + // Return energy
  500 + return energy;
  501 +
  502 +}
  503 +#else
391 504 GEnergy GCTAEdisp2D::mc(GRan& ran,
392 505 const double& logEsrc,
393 506 const double& theta,
... ... @@ -440,8 +553,6 @@ GEnergy GCTAEdisp2D::mc(GRan&amp; ran,
440 553 index = low;
441 554 }
442 555  
443   - //std::cout << "p = " << p << std::endl;
444   - //std::cout << "index = " << index << std::endl;
445 556  
446 557 // Interpolate EobsOverEsrc value
447 558 double EobsOverEsrc = 0.0;
... ... @@ -459,12 +570,11 @@ GEnergy GCTAEdisp2D::mc(GRan&amp; ran,
459 570 // Set energy
460 571 GEnergy energy;
461 572 energy.TeV(EobsOverEsrc * Esrc);
462   - //std::cout << "eobs = " << EobsOverEsrc * Esrc << std::endl;
463 573  
464 574 // Return energy
465 575 return energy;
466 576 }
467   -
  577 +#endif
468 578  
469 579 /***********************************************************************//**
470 580 * @brief Return observed energy interval that contains the energy dispersion.
... ... @@ -486,57 +596,42 @@ GEbounds GCTAEdisp2D::ebounds_obs(const double&amp; logEsrc,
486 596 const double& zenith,
487 597 const double& azimuth) const
488 598 {
489   - // Set epsilon
490   - const double eps = 1.0e-06;
  599 + // Compute only if parameters changed
  600 + if (!m_ebounds_obs_computed || theta != m_theta) {
491 601  
492   - // Initialise results
493   - double logEobsMin = -10.0;
494   - double logEobsMax = 30.0;
495   - bool minFound = false;
496   - bool maxFound = false;
  602 + m_ebounds_obs_computed = true;
  603 + m_theta = theta;
497 604  
498   - // Find boundaries
499   - for (int i = 0; i < m_edisp.axis(1); ++i) {
  605 + // Compute ebounds_obs
  606 + compute_ebounds_obs(theta, phi, zenith, azimuth);
500 607  
501   - // Compute value EobsOverEsrc
502   - double EobsOverEsrc = 0.5 * (m_edisp.axis_hi(1, i) +
503   - m_edisp.axis_lo(1, i));
  608 + }
504 609  
505   - // Find first non-negligible matrix term
506   - if (!minFound && m_edisp(0, logEsrc, EobsOverEsrc, theta) >= eps) {
507   - minFound = true;
508   - logEobsMin = std::log10(EobsOverEsrc) + logEsrc;
509   - }
  610 + // Search index only if logEsrc has changed
  611 + if (logEsrc != m_logEsrc) {
510 612  
511   - // Find last non-negligible matrix term
512   - else if (minFound && !maxFound &&
513   - m_edisp(0, logEsrc, EobsOverEsrc, theta) < eps) {
514   - maxFound = true;
515   - logEobsMax = std::log10(EobsOverEsrc) + logEsrc;
516   - }
  613 + m_logEsrc = logEsrc;
517 614  
518   - // Continue searching
519   - else if (minFound && maxFound &&
520   - m_edisp(0, logEsrc, EobsOverEsrc, theta) >= eps) {
521   - maxFound = false;
  615 + // Find right index with bisection
  616 + int low = 0;
  617 + int high = m_ebounds_obs.size();
  618 + while ((high-low) > 1) {
  619 + int mid = (low+high) / 2;
  620 + double e = m_edisp.axis_lo(0, mid);
  621 + if (logEsrc < std::log10(e)) {
  622 + high = mid;
  623 + }
  624 + else {
  625 + low = mid;
  626 + }
522 627 }
  628 + // Index found
  629 + m_index_obs = low;
523 630  
524   - } // endfor: looped over EobsOverEsrc
525   -
526   - // If energy dispersion has never become negligible until end of loop,
527   - // reset logEobsMax
528   - if (!maxFound) {
529   - logEobsMax = 30.0;
530 631 }
531 632  
532   - // Compute energy boundaries
533   - GEnergy emin;
534   - GEnergy emax;
535   - emin.log10TeV(logEobsMin);
536   - emax.log10TeV(logEobsMax);
537   -
538 633 // Return energy boundaries
539   - return (GEbounds(emin, emax));
  634 + return m_ebounds_obs[m_index_obs];
540 635 }
541 636  
542 637  
... ... @@ -560,58 +655,48 @@ GEbounds GCTAEdisp2D::ebounds_src(const double&amp; logEobs,
560 655 const double& zenith,
561 656 const double& azimuth) const
562 657 {
563   - // Set epsilon
564   - const double eps = 1.0e-06;
  658 + // Compute only if parameters changed
  659 + if (!m_ebounds_src_computed || theta != m_theta) {
565 660  
566   - // Initialise results
567   - double logEsrcMin = -10.0;
568   - double logEsrcMax = 30.0;
569   - bool minFound = false;
570   - bool maxFound = false;
  661 + m_ebounds_src_computed = true;
  662 + m_theta = theta;
571 663  
572   - // Find boundaries
573   - for (int i = 0; i < m_edisp.axis(0); ++i) {
  664 + // Compute ebounds_src
  665 + compute_ebounds_src(theta, phi, zenith, azimuth);
574 666  
575   - // Compute value of corresponding logEsrc
576   - double logEsrc = 0.5 * (m_edisp.axis_hi(0, i) +
577   - m_edisp.axis_lo(0, i));
578   - double EobsOverEsrc = std::exp((logEobs-logEsrc)*gammalib::ln10);
  667 + }
579 668  
580   - // Find first non-negligible matrix term
581   - if (!minFound && m_edisp(0, logEsrc, EobsOverEsrc, theta) >= eps) {
582   - minFound = true;
583   - logEsrcMin = logEsrc;
584   - }
  669 + // Search index only if logEobs has changed
  670 + if (logEobs != m_logEobs) {
585 671  
586   - // Find last non-negligible matrix term
587   - else if (minFound && !maxFound &&
588   - m_edisp(0, logEsrc, EobsOverEsrc, theta) < eps) {
589   - maxFound = true;
590   - logEsrcMax = logEsrc;
591   - }
  672 + m_logEobs = logEobs;
592 673  
593   - // Continue
594   - else if (minFound && maxFound &&
595   - m_edisp(0, logEsrc, EobsOverEsrc, theta) >= eps) {
596   - maxFound = false;
597   - }
  674 + // Find right index with bisection
  675 + int low = 0;
  676 + int high = m_ebounds_src.size();
598 677  
599   - } // endfor: looped over true energy
  678 + while ((high-low) > 1) {
  679 + int mid = (low+high) / 2;
600 680  
601   - // If energy dispersion has never become negligible until end of loop,
602   - // reset logEobsMax
603   - if (!maxFound) {
604   - logEsrcMax = 30.0;
605   - }
  681 + double e = m_edisp.axis_lo(0, mid);
  682 + //double e = m_eobs_axis[mid];
  683 + if (logEobs < std::log10(e)) {
  684 + high = mid;
  685 + }
  686 + else {
  687 + low = mid;
  688 + }
  689 + }
  690 + // Index found
  691 + m_index_src = low;
  692 + //m_index_src = low * m_edisp.axis(1) / m_edisp.axis(0);
  693 + //std::cout << "emin = " << m_ebounds_src[m_index_src].emin().TeV() << std::endl;
  694 + //std::cout << "emax = " << m_ebounds_src[m_index_src].emax().TeV() << std::endl;
606 695  
607   - // Compute energy boundaries
608   - GEnergy emin;
609   - GEnergy emax;
610   - emin.log10TeV(logEsrcMin);
611   - emax.log10TeV(logEsrcMax);
  696 + }
612 697  
613 698 // Return energy boundaries
614   - return (GEbounds(emin, emax));
  699 + return m_ebounds_src[m_index_src];
615 700 }
616 701  
617 702  
... ... @@ -669,15 +754,30 @@ std::string GCTAEdisp2D::print(const GChatter&amp; chatter) const
669 754 ***************************************************************************/
670 755 void GCTAEdisp2D::init_members(void)
671 756 {
  757 +
672 758 // Initialise members
673 759 m_filename.clear();
674 760 m_edisp.clear();
675 761  
676 762 // Initialise Monte Carlo cache
677 763 m_theta = 0.0;
  764 + m_logEsrc =-30.0;
  765 + m_logEobs =-30.0;
  766 + m_fmax = -1.0;
  767 + m_index_obs = 0;
  768 + m_index_src = 0;
678 769 m_cdf_computed = false;
  770 + m_ebounds_obs_computed = false;
  771 + m_ebounds_src_computed = false;
  772 +
679 773 m_cumul.clear();
680 774  
  775 + // Initialize ebounds
  776 + m_ebounds_obs.clear();
  777 + m_ebounds_src.clear();
  778 +
  779 + m_eobs_axis.clear();
  780 +
681 781 // Return
682 782 return;
683 783 }
... ... @@ -691,8 +791,19 @@ void GCTAEdisp2D::init_members(void)
691 791 void GCTAEdisp2D::copy_members(const GCTAEdisp2D& edisp)
692 792 {
693 793 // Copy members
694   - m_filename = edisp.m_filename;
695   - m_edisp = edisp.m_edisp;
  794 + m_filename = edisp.m_filename;
  795 + m_edisp = edisp.m_edisp;
  796 + m_theta = edisp.m_theta;
  797 + m_logEsrc = edisp.m_logEsrc;
  798 + m_logEobs = edisp.m_logEobs;
  799 + m_index_obs = edisp.m_index_obs;
  800 + m_index_src = edisp.m_index_src;
  801 +
  802 + m_ebounds_obs = edisp.m_ebounds_obs;
  803 + m_ebounds_src = edisp.m_ebounds_src;
  804 + m_cumul = edisp.m_cumul;
  805 + m_eobs_axis = edisp.m_eobs_axis;
  806 + m_temp = edisp.m_temp;
696 807  
697 808 // Return
698 809 return;
... ... @@ -709,9 +820,239 @@ void GCTAEdisp2D::free_members(void)
709 820 }
710 821  
711 822 /***********************************************************************//**
  823 + * @brief Compute ebounds_obs vector
  824 + *
  825 + ***************************************************************************/
  826 +void GCTAEdisp2D::compute_ebounds_obs(const double& theta,
  827 + const double& phi,
  828 + const double& zenith,
  829 + const double& azimuth) const
  830 +{
  831 + // Clear ebounds_obs vector
  832 + m_ebounds_obs.clear();
  833 +
  834 + // Set epsilon
  835 + const double eps = 1.0e-06;
  836 +
  837 + // Loop over Esrc
  838 + for (int isrc = 0; isrc < m_edisp.axis(0); ++isrc) {
  839 +
  840 + // Set Esrc
  841 + double Esrc = 0.5 * (m_edisp.axis_hi(0, isrc) +
  842 + m_edisp.axis_lo(0, isrc));
  843 + double logEsrc = std::log10(Esrc);
  844 +
  845 + // Initialise results
  846 + double logEobsMin = -30.0;
  847 + double logEobsMax = 10.0;
  848 + bool minFound = false;
  849 + bool maxFound = false;
  850 +
  851 + // Find boundaries, loop over EobsOverEsrc
  852 + for (int i = 0; i < m_edisp.axis(1); ++i) {
  853 +
  854 + // Compute value EobsOverEsrc
  855 + double EobsOverEsrc = 0.5 * (m_edisp.axis_hi(1, i) +
  856 + m_edisp.axis_lo(1, i));
  857 +
  858 + // Find first non-negligible matrix term
  859 + if (!minFound && m_edisp(0, logEsrc, EobsOverEsrc, theta) >= eps) {
  860 + minFound = true;
  861 + logEobsMin = std::log10(EobsOverEsrc * Esrc);
  862 + }
  863 +
  864 + // Find last non-negligible matrix term
  865 + else if (minFound && !maxFound &&
  866 + m_edisp(0, logEsrc, EobsOverEsrc, theta) < eps) {
  867 + maxFound = true;
  868 + logEobsMax = std::log10(EobsOverEsrc * Esrc);
  869 + }
  870 +
  871 + // Continue searching
  872 + else if (minFound && maxFound &&
  873 + m_edisp(0, logEsrc, EobsOverEsrc, theta) >= eps) {
  874 + maxFound = false;
  875 + }
  876 +
  877 + } // endfor: looped over EobsOverEsrc
  878 +
  879 +
  880 + // If energy dispersion has never become negligible until end of loop,
  881 + // reset logEobsMax
  882 + if (!maxFound) {
  883 + logEobsMax = 10.0;
  884 + }
  885 +
  886 + // Set energy boundaries
  887 + GEnergy emin;
  888 + GEnergy emax;
  889 + emin.log10TeV(logEobsMin);
  890 + emax.log10TeV(logEobsMax);
  891 +
  892 + // Add energy boundaries to vector
  893 + m_ebounds_obs.push_back(GEbounds(emin, emax));
  894 +
  895 + }
  896 +
  897 +
  898 + // Return
  899 + return;
  900 +}
  901 +
  902 +/***********************************************************************//**
  903 + * @brief Compute ebounds_src vector
  904 + *
  905 + ***************************************************************************/
  906 +void GCTAEdisp2D::compute_ebounds_src(const double& theta,
  907 + const double& phi,
  908 + const double& zenith,
  909 + const double& azimuth) const
  910 +{
  911 +
  912 + // Clear ebounds_src vector
  913 + m_ebounds_src.clear();
  914 +
  915 + // Set epsilon
  916 + const double eps = 1.0e-06;
  917 +
  918 +/*
  919 + // Loop over EobsOverEsrc
  920 + for (int iobs = 0; iobs < m_edisp.axis(1); ++iobs) {
  921 +
  922 + // Set EobsOverEsrc
  923 + double EobsOverEsrc = 0.5 * (m_edisp.axis_hi(1, iobs) +
  924 + m_edisp.axis_lo(1, iobs));
  925 +
  926 + // Initialise results
  927 + double logEsrcMin = -10.0;
  928 + double logEsrcMax = 30.0;
  929 + bool minFound = false;
  930 + bool maxFound = false;
  931 +
  932 + // Find boundaries, loop over Esrc
  933 + for (int i = 0; i < m_edisp.axis(0); ++i) {
  934 +
  935 + // Compute value of corresponding logEsrc
  936 + double Esrc = 0.5 * (m_edisp.axis_hi(0, i) +
  937 + m_edisp.axis_lo(0, i));
  938 + double logEsrc = std::log10(Esrc);
  939 +
  940 + // Find first non-negligible matrix term
  941 + if (!minFound && m_edisp(0, logEsrc, EobsOverEsrc, theta) >= eps) {
  942 + minFound = true;
  943 + logEsrcMin = logEsrc;
  944 + }
  945 +
  946 + // Find last non-negligible matrix term
  947 + else if (minFound && !maxFound &&
  948 + m_edisp(0, logEsrc, EobsOverEsrc, theta) < eps) {
  949 + maxFound = true;
  950 + logEsrcMax = logEsrc;
  951 + }
  952 +
  953 + // Continue
  954 + else if (minFound && maxFound &&
  955 + m_edisp(0, logEsrc, EobsOverEsrc, theta) >= eps) {
  956 + maxFound = false;
  957 + }
  958 +
  959 + } // endfor: looped over true energy
  960 +
  961 + // If energy dispersion has never become negligible until end of loop,
  962 + // reset logEobsMax
  963 + if (!maxFound) {
  964 + logEsrcMax = 30.0;
  965 + }
  966 +
  967 + // Set energy boundaries
  968 + GEnergy emin;
  969 + GEnergy emax;
  970 + emin.log10TeV(logEsrcMin);
  971 + emax.log10TeV(logEsrcMax);
  972 +
  973 + // Add energy boundaries to vector
  974 + m_ebounds_src.push_back(GEbounds(emin, emax));
  975 +
  976 + }
  977 +
  978 + const int step = m_edisp.axis(1) / m_edisp.axis(0);
  979 +
  980 + // Loop over EobsOverEsrc to create Eobs axis
  981 + for (int i = 0; i < m_edisp.axis(0); ++i) {
  982 + m_eobs_axis.push_back(m_edisp.axis_lo(1, i*step) * m_edisp.axis_lo(0, i));
  983 + }
  984 +*/
  985 +
  986 + // Loop over Eobs
  987 + for (int iobs = 0; iobs < m_edisp.axis(0); ++iobs) {
  988 +
  989 + double Eobs = 0.5 * (m_edisp.axis_hi(0, iobs) +
  990 + m_edisp.axis_lo(0, iobs));
  991 + double logEobs = std::log10(Eobs);
  992 +
  993 + // Initialise results
  994 + double logEsrcMin = -10.0;
  995 + double logEsrcMax = 30.0;
  996 + bool minFound = false;
  997 + bool maxFound = false;
  998 +
  999 + // Find boundaries, loop over Esrc
  1000 + for (int isrc = 0; isrc < m_edisp.axis(0); ++isrc) {
  1001 +
  1002 + double Esrc = 0.5 * (m_edisp.axis_hi(0, isrc) +
  1003 + m_edisp.axis_lo(0, isrc));
  1004 + double logEsrc = std::log10(Esrc);
  1005 +
  1006 + // Find first non-negligible matrix term
  1007 + if (!minFound && operator()(logEobs, logEsrc, theta) >= eps) {
  1008 + minFound = true;
  1009 + logEsrcMin = logEsrc;
  1010 + }
  1011 +
  1012 + // Find last non-negligible matrix term
  1013 + else if (minFound && !maxFound &&
  1014 + operator()(logEobs, logEsrc, theta) < eps) {
  1015 + maxFound = true;
  1016 + logEsrcMax = logEsrc;
  1017 + }
  1018 +
  1019 + // Continue
  1020 + else if (minFound && maxFound &&
  1021 + operator()(logEobs, logEsrc, theta) >= eps) {
  1022 + maxFound = false;
  1023 + }
  1024 +
  1025 + } // endfor: looped over true energy
  1026 +
  1027 + // If energy dispersion has never become negligible until end of loop,
  1028 + // reset logEobsMax
  1029 + if (!maxFound) {
  1030 + logEsrcMax = 30.0;
  1031 + }
  1032 +
  1033 +
  1034 + // Set energy boundaries
  1035 + GEnergy emin;
  1036 + GEnergy emax;
  1037 + emin.log10TeV(logEsrcMin);
  1038 + emax.log10TeV(logEsrcMax);
  1039 +
  1040 + // Add energy boundaries to vector
  1041 + m_ebounds_src.push_back(GEbounds(emin, emax));
  1042 +
  1043 + } // endfor : loopend over observed energy
  1044 +
  1045 +
  1046 + // Return
  1047 + return;
  1048 +}
  1049 +
  1050 +
  1051 +/***********************************************************************//**
712 1052 * @brief Compute cumulative probability
713 1053 *
714 1054 ***************************************************************************/
  1055 +#if defined(G_CDF_METHOD)
715 1056 void GCTAEdisp2D::compute_cumul(const double& theta,
716 1057 const double& phi,
717 1058 const double& zenith,
... ... @@ -737,39 +1078,10 @@ void GCTAEdisp2D::compute_cumul(const double&amp; theta,
737 1078 // Initialize cumulative probability
738 1079 double sum = 0.0;
739 1080  
740   - // Loop over Eobs from 0.1 TeV to 100.0 TeV
741   - //int nsteps = (int) (m_edisp.axis_hi(1, m_edisp.axis(1)-1) * Esrc / deltaEobs);
742   - int nsteps = 10000;
743   - double deltaEobs = 100.0/nsteps;
744   - double Eobs = 0.1;
745   -
746   - for (int iobs = 0; iobs < nsteps; ++iobs) {
747   -
748   - double EobsOverEsrc = Eobs/Esrc;
749   -
750   - // Compute cumulative probability
751   - double add = m_edisp(0, std::log10(Esrc), EobsOverEsrc, theta) *
752   - deltaEobs / Eobs;
753   -
754   - // Add to sum (and keep sum lower or equal to 1)
755   - sum = sum+add >= 1.0 ? 1.0 : sum+add;
756   -
757   - // Create pair containing EobsOverEsrc and cumulated probability
758   - std::pair<double, double> pair(EobsOverEsrc, sum);
759   -
760   - // Add to vector
761   - m_cumul[isrc].push_back(pair);
762   -
763   - // Compute new Eobs value
764   - Eobs += deltaEobs;
765   -
766   - } // endfor: looped over Eobs
767   -
768   - /* // Loop over EobsOverEtrue
769   -
770 1081 // Divide bin into n sub-bins
771 1082 const int n = 10;
772 1083  
  1084 + // Loop over EobsOverEtrue
773 1085 for (int imigra = 0; imigra < m_edisp.axis(1); ++imigra) {
774 1086  
775 1087 // Compute delta(Eobs/Etrue) and Eobs/Etrue values
... ... @@ -791,16 +1103,13 @@ void GCTAEdisp2D::compute_cumul(const double&amp; theta,
791 1103 // Create pair containing EobsOverEsrc and cumulated probability
792 1104 std::pair<double, double> pair(EobsOverEsrc, sum);
793 1105  
794   - //if(Esrc == 10.0 && EobsOverEsrc < 2.0)
795   - //std::cout << "Eobs/Esrc = " << EobsOverEsrc << ", add = " << add << ", sum = " << sum << std::endl;
796   -
797 1106 // Add to vector
798 1107 m_cumul[isrc].push_back(pair);
799 1108  
800 1109 }
801 1110  
802 1111 } // endfor: looped over EobsOverEsrc
803   - */
  1112 +
804 1113 }
805 1114  
806 1115 }
... ... @@ -808,3 +1117,4 @@ void GCTAEdisp2D::compute_cumul(const double&amp; theta,
808 1117 // Return
809 1118 return;
810 1119 }
  1120 +#endif
... ...
inst/cta/src/GCTAEdispPerfTable.cpp
  View file @cee203a
... ... @@ -245,7 +245,7 @@ GCTAEdispPerfTable* GCTAEdispPerfTable::clone(void) const
245 245 * performance table. The energy resolution is stored in the 5th column
246 246 * of the performance table as RMS(ln(Eest/Etrue)). The method converts
247 247 * this internally to a sigma value by multiplying the stored values by
248   - * 0.434294481903.
  248 + * 1/ln(10).
249 249 ***************************************************************************/
250 250 void GCTAEdispPerfTable::load(const std::string& filename)
251 251 {
... ... @@ -290,8 +290,15 @@ void GCTAEdispPerfTable::load(const std::string&amp; filename)
290 290  
291 291 // Push elements in node array and vector
292 292 m_logE.append(gammalib::todouble(elements[0]));
  293 + // The energy resolution is stored as RMS(ln(Eest/Etrue))
293 294 m_sigma.push_back(gammalib::todouble(elements[4]) * gammalib::inv_ln10);
294 295  
  296 + // The energy resolution is stored as RMS(Eest/Etrue)
  297 + //m_sigma.push_back(std::log10(gammalib::todouble(elements[4])));
  298 +
  299 + // The energy resolution is stored as RMS(log10(Eest/Etrue))
  300 + //m_sigma.push_back(gammalib::todouble(elements[4]));
  301 +
295 302 } // endwhile: looped over lines
296 303  
297 304 // Close file
... ... @@ -309,7 +316,7 @@ void GCTAEdispPerfTable::load(const std::string&amp; filename)
309 316 * @brief Simulate energy dispersion
310 317 *
311 318 * @param[in] ran Random number generator.
312   - * @param[in] logE Log10 of the true photon energy (TeV).
  319 + * @param[in] logEsrc Log10 of the true photon energy (TeV).
313 320 * @param[in] theta Offset angle in camera system (rad). Not used.
314 321 * @param[in] phi Azimuth angle in camera system (rad). Not used.
315 322 * @param[in] zenith Zenith angle in Earth system (rad). Not used.
... ... @@ -319,17 +326,17 @@ void GCTAEdispPerfTable::load(const std::string&amp; filename)
319 326 * m_par_sigma around @p logE.
320 327 ***************************************************************************/
321 328 GEnergy GCTAEdispPerfTable::mc(GRan& ran,
322   - const double& logE,
  329 + const double& logEsrc,
323 330 const double& theta,
324 331 const double& phi,
325 332 const double& zenith,
326 333 const double& azimuth) const
327 334 {
328 335 // Update the parameter cache
329   - update(logE);
  336 + update(logEsrc);
330 337  
331 338 // Draw log observed energy in TeV
332   - double logEobs = m_par_sigma * ran.normal() + logE;
  339 + double logEobs = m_par_sigma * ran.normal() + logEsrc;
333 340  
334 341 // Set energy
335 342 GEnergy energy;
... ... @@ -530,7 +537,7 @@ void GCTAEdispPerfTable::update(const double&amp; logE) const
530 537  
531 538 // Save energy
532 539 m_par_logE = logE;
533   -
  540 +
534 541 // Determine Gaussian sigma and pre-compute Gaussian parameters
535 542 m_par_sigma = m_logE.interpolate(logE, m_sigma);
536 543 m_par_scale = gammalib::inv_sqrt2pi / m_par_sigma;
... ...
inst/cta/src/GCTAEdispRmf.cpp
  View file @cee203a
... ... @@ -42,6 +42,8 @@
42 42 /* __ Coding definitions _________________________________________________ */
43 43 //#define G_OLD_MC_CODE
44 44 //#define G_MC_REJECTION
  45 +#define G_NEW_REJECTION
  46 +//#define G_CDF_METHOD
45 47  
46 48 /* __ Debug definitions __________________________________________________ */
47 49  
... ... @@ -271,6 +273,40 @@ void GCTAEdispRmf::load(const std::string&amp; filename)
271 273 * more efficient. Either implement a rejection method for smoothness
272 274 * or introduce subsampling for CDF cache.
273 275 ***************************************************************************/
  276 +#if defined(G_NEW_REJECTION)
  277 +GEnergy GCTAEdispRmf::mc(GRan& ran,
  278 + const double& logEsrc,
  279 + const double& theta,
  280 + const double& phi,
  281 + const double& zenith,
  282 + const double& azimuth) const
  283 +{
  284 +
  285 + // Set true energy
  286 + GEnergy energy;
  287 + energy.log10TeV(logEsrc);
  288 +
  289 + // Get boundaries for rejection method
  290 + double emin = ebounds_obs(logEsrc, theta, phi, zenith, azimuth).emin().log10TeV();
  291 + double emax = ebounds_obs(logEsrc, theta, phi, zenith, azimuth).emax().log10TeV();
  292 + double fmax = m_matrix(m_rmf.itruemax(), m_rmf.imeasmax());
  293 +
  294 + // Find energy by rejection method
  295 + double ewidth = emax - emin;
  296 + double e = emin;
  297 + double f = 0.0;
  298 + double ftest = 1.0;
  299 + while (ftest > f) {
  300 + e = emin + ewidth * ran.uniform();
  301 + f = operator()(e, logEsrc, theta, phi, zenith, azimuth);
  302 + ftest = ran.uniform() * fmax;
  303 + }
  304 + energy.log10TeV(e);
  305 +
  306 + // Return energy
  307 + return energy;
  308 +}
  309 +#endif
274 310 #if defined(G_OLD_MC_CODE)
275 311 GEnergy GCTAEdispRmf::mc(GRan& ran,
276 312 const double& logEsrc,
... ... @@ -356,7 +392,8 @@ GEnergy GCTAEdispRmf::mc(GRan&amp; ran,
356 392 // Return energy
357 393 return energy;
358 394 }
359   -#else
  395 +#endif
  396 +#if defined(G_CDF_METHOD)
360 397 GEnergy GCTAEdispRmf::mc(GRan& ran,
361 398 const double& logEsrc,
362 399 const double& theta,
... ... @@ -449,14 +486,73 @@ GEbounds GCTAEdispRmf::ebounds_obs(const double&amp; logEsrc,
449 486 const double& zenith,
450 487 const double& azimuth) const
451 488 {
452   - // Set true energy
453   - GEnergy etrue;
454   - etrue.log10TeV(logEsrc);
455 489  
456   - // Return measured energy boundaries
457   - return (m_rmf.emeasured(etrue));
  490 + // Compute only if parameters changed
  491 + if (!m_ebounds_obs_computed || theta != m_theta) {
  492 +
  493 + m_ebounds_obs_computed = true;
  494 + m_theta = theta;
  495 +
  496 + // Compute ebounds_obs
  497 + compute_ebounds_obs(theta, phi, zenith, azimuth);
  498 +
  499 + }
  500 +
  501 + // Search index only if logEsrc has changed
  502 + if (logEsrc != m_logEsrc) {
  503 +
  504 + m_logEsrc = logEsrc;
  505 +
  506 + // Find right index with bisection
  507 + int low = 0;
  508 + int high = m_ebounds_obs.size();
  509 + while ((high-low) > 1) {
  510 + int mid = (low+high) / 2;
  511 + double e = m_rmf.etrue().emin(mid).log10TeV();
  512 + if (logEsrc < e) {
  513 + high = mid;
  514 + }
  515 + else {
  516 + low = mid;
  517 + }
  518 + }
  519 + // Index found
  520 + m_index_obs = low;
  521 +
  522 + }
  523 +
  524 + // Return energy boundaries
  525 + return m_ebounds_obs[m_index_obs];
  526 +
458 527 }
459 528  
  529 +/***********************************************************************//**
  530 + * @brief Compute ebounds_obs vector
  531 + *
  532 + ***************************************************************************/
  533 +void GCTAEdispRmf::compute_ebounds_obs(const double& theta,
  534 + const double& phi,
  535 + const double& zenith,
  536 + const double& azimuth) const
  537 +{
  538 +
  539 + // Loop over Etrue
  540 + for (int itrue = 0; itrue < m_rmf.ntrue(); ++itrue) {
  541 +
  542 + double logEsrc = m_rmf.etrue().emin(itrue).log10TeV();
  543 +
  544 + // Set true energy
  545 + GEnergy etrue;
  546 + etrue.log10TeV(logEsrc);
  547 +
  548 + // Add ebounds to vector
  549 + m_ebounds_obs.push_back(m_rmf.emeasured(etrue));
  550 + }
  551 +
  552 + // Return
  553 + return;
  554 +
  555 +}
460 556  
461 557 /***********************************************************************//**
462 558 * @brief Return true energy interval that contains the energy dispersion.
... ... @@ -476,14 +572,75 @@ GEbounds GCTAEdispRmf::ebounds_src(const double&amp; logEobs,
476 572 const double& zenith,
477 573 const double& azimuth) const
478 574 {
479   - // Set measured energy
480   - GEnergy emeasured;
481   - emeasured.log10TeV(logEobs);
482 575  
483   - // Return true energy boundaries
484   - return (m_rmf.etrue(emeasured));
  576 + // Compute only if parameters changed
  577 + if (!m_ebounds_src_computed || theta != m_theta) {
  578 +
  579 + m_ebounds_src_computed = true;
  580 + m_theta = theta;
  581 +
  582 + // Compute ebounds_obs
  583 + compute_ebounds_src(theta, phi, zenith, azimuth);
  584 +
  585 + }
  586 +
  587 + // Search index only if logEobs has changed
  588 + if (logEobs != m_logEobs) {
  589 +
  590 + m_logEobs = logEobs;
  591 +
  592 + // Find right index with bisection
  593 + int low = 0;
  594 + int high = m_ebounds_src.size();
  595 + while ((high-low) > 1) {
  596 + int mid = (low+high) / 2;
  597 + double e = m_rmf.emeasured().emin(mid).log10TeV();
  598 + if (logEobs < e) {
  599 + high = mid;
  600 + }
  601 + else {
  602 + low = mid;
  603 + }
  604 + }
  605 + // Index found
  606 + m_index_src = low;
  607 + //std::cout << "emin = " << m_ebounds_src[m_index_src].emin().TeV() << std::endl;
  608 + //std::cout << "emax = " << m_ebounds_src[m_index_src].emax().TeV() << std::endl;
  609 +
  610 + }
  611 +
  612 + // Return energy boundaries
  613 + return m_ebounds_src[m_index_src];
  614 +
485 615 }
486 616  
  617 +/***********************************************************************//**
  618 + * @brief Compute ebounds_src vector
  619 + *
  620 + ***************************************************************************/
  621 +void GCTAEdispRmf::compute_ebounds_src(const double& theta,
  622 + const double& phi,
  623 + const double& zenith,
  624 + const double& azimuth) const
  625 +{
  626 +
  627 + // Loop over Eobs
  628 + for (int imeas = 0; imeas < m_rmf.nmeasured(); ++imeas) {
  629 +
  630 + double logEobs = m_rmf.emeasured().emin(imeas).log10TeV();
  631 +
  632 + // Set measured energy
  633 + GEnergy emeasured;
  634 + emeasured.log10TeV(logEobs);
  635 +
  636 + // Add ebounds to vector
  637 + m_ebounds_src.push_back(m_rmf.emeasured(emeasured));
  638 + }
  639 +
  640 + // Return
  641 + return;
  642 +
  643 +}
487 644  
488 645 /***********************************************************************//**
489 646 * @brief Print RMF information
... ... @@ -560,8 +717,17 @@ void GCTAEdispRmf::init_members(void)
560 717 m_mc_measured_cdf.clear();
561 718 m_cdf_computed = false;
562 719 m_theta = 0.0;
  720 + m_logEsrc = -30.0;
  721 + m_logEobs = -30.0;
563 722 m_cumul.clear();
564 723  
  724 + m_ebounds_obs_computed = false;
  725 + m_ebounds_obs.clear();
  726 + m_index_obs = 0;
  727 + m_ebounds_src_computed = false;
  728 + m_ebounds_src.clear();
  729 + m_index_src = 0;
  730 +
565 731 // Return
566 732 return;
567 733 }
... ...
inst/cta/src/GCTAResponseTable.cpp
  View file @cee203a
... ... @@ -1753,6 +1753,7 @@ void GCTAResponseTable::update(const double&amp; arg1, const double&amp; arg2) const
1753 1753 void GCTAResponseTable::update(const double& arg1, const double& arg2,
1754 1754 const double& arg3) const
1755 1755 {
  1756 +
1756 1757 // Get pointers to node arrays (circumvent const correctness)
1757 1758 const GNodeArray* nodes1 = &(m_axis_nodes[0]);
1758 1759 const GNodeArray* nodes2 = &(m_axis_nodes[1]);
... ...
inst/cta/test/caldb/edisp_matrix.fits
View file @cee203a
No preview for this file type
inst/cta/test/example_sim_edisp.py
  View file @cee203a
... ... @@ -79,17 +79,21 @@ def sim_edisp(edisp, etrue, eobs_max=40.0, ebins=1000, nmc=100000):
79 79 sum = 0.0
80 80 exp_edisp = []
81 81 cumul = []
  82 + #fmax = 0.0
82 83 for i in range(ebins):
83 84 eobs = eobs_axis[i]
84 85 if eobs <= edisp.ebounds_obs(logEtrue).emax().TeV():
85 86 dLogEobs = math.log10(eobs_axis[i]) - math.log10(eobs_axis[i-1])
86 87 value = edisp(math.log10(eobs), logEtrue) * dLogEobs * nmc
  88 + #if fmax < value / dLogEobs / nmc:
  89 + # fmax = value / dLogEobs / nmc
87 90 else:
88 91 value = 0.0
89 92 sum += value
90 93 exp_edisp.append(value)
91 94 cumul.append(sum/nmc)
92 95 print(sum/nmc)
  96 + #print(fmax)
93 97  
94 98 # Create figure
95 99 fig, ax1 = plt.subplots()
... ... @@ -105,13 +109,13 @@ def sim_edisp(edisp, etrue, eobs_max=40.0, ebins=1000, nmc=100000):
105 109 ax1.set_xlabel("Observed energy (TeV)")
106 110 ax1.set_ylabel("Number of counts")
107 111  
108   - ax2 = ax1.twinx()
  112 + #ax2 = ax1.twinx()
109 113  
110 114 # Plot cumulative dispersion
111   - ax2.plot(eobs_axis, cumul, 'm')
  115 + #ax2.plot(eobs_axis, cumul, 'm')
112 116  
113 117 # Set axis
114   - ax2.set_ylabel("Cumulative probability")
  118 + #ax2.set_ylabel("Cumulative probability")
115 119  
116 120 # Notify
117 121 print("PLEASE CLOSE WINDOW TO CONTINUE ...")
... ... @@ -140,7 +144,7 @@ if __name__ == &#39;__main__&#39;:
140 144 # Load edisp
141 145 #edisp = gammalib.GCTAEdispRmf("./caldb/dc1/rmf.fits")
142 146 #edisp = gammalib.GCTAEdispPerfTable("./caldb/cta_dummy_irf.dat")
143   - edisp = gammalib.GCTAEdisp2D("./caldb/edisp_matrix.fits")
  147 + edisp = gammalib.GCTAEdisp2D("/home/forest/git/gammalib/inst/cta/test/caldb/edisp_matrix.fits")
144 148  
145 149 # Simulate Edisp
146 150 sim_edisp(edisp, 20.0) # Energy in TeV
... ...
inst/cta/test/test_CTA.cpp
  View file @cee203a
... ... @@ -611,12 +611,6 @@ void TestGCTAResponse::test_response_edisp_RMF(void)
611 611 // Test normalisation
612 612 test_edisp_integration(edisp, 0.5, 10.0);
613 613  
614   - std::cout << "Eobs = 1 TeV : " << edisp(0.0, 0.0) << std::endl;
615   - std::cout << "Eobs = 10 TeV : " << edisp(1.0, 0.0) << std::endl;
616   -
617   - std::cout << "Eobs = 10 TeV : " << edisp(1.0, 1.0) << std::endl;
618   - std::cout << "Eobs = 100 TeV : " << edisp(2.0, 1.0) << std::endl;
619   -
620 614 // Return
621 615 return;
622 616 }
... ... @@ -650,23 +644,8 @@ void TestGCTAResponse::test_response_edisp_2D(void)
650 644 GCTAEdisp2D edisp(cta_edisp_2D);
651 645  
652 646 // Test normalisation
653   - test_edisp_integration(edisp);
  647 + test_edisp_integration(edisp, 10.0, 50.0);
654 648  
655   - // Test ebounds_obs and ebounds_src methods
656   - test_try("GCTAEdisp2D ebounds_src and ebounds_obs methods");
657   - try {
658   - GCTAEdisp2D edisp(cta_edisp_2D);
659   - test_try_success();
660   - }
661   - catch (std::exception &e) {
662   - test_try_failure(e);
663   - }
664   -
665   - std::cout << "Eobs = 1 TeV : " << edisp(0.0, 0.0) << std::endl;
666   - std::cout << "Eobs = 10 TeV : " << edisp(1.0, 0.0) << std::endl;
667   -
668   - std::cout << "Eobs = 10 TeV : " << edisp(1.0, 1.0) << std::endl;
669   - std::cout << "Eobs = 100 TeV : " << edisp(2.0, 1.0) << std::endl;
670 649  
671 650 // Return
672 651 return;
... ... @@ -997,14 +976,14 @@ void TestGCTAResponse::test_edisp_integration(const GCTAEdisp&amp; edisp,
997 976 const double& e_src_max)
998 977 {
999 978  
1000   - // Loop over source energies
1001   - for (double e_src = e_src_min; e_src < e_src_max; e_src *= 2.0) {
  979 + // Loop over source energies
  980 + for (double e_src = e_src_min; e_src < e_src_max; e_src *= 2.0) {
1002 981  
1003   - // Compute log10 of true energy
1004   - double log10_e_src = std::log10(e_src);
  982 + // Compute log10 of true energy
  983 + double log10_e_src = std::log10(e_src);
1005 984  
1006   - // Retrieve boundaries in observed energy
1007   - GEbounds ebounds = edisp.ebounds_obs(log10_e_src);
  985 + // Retrieve boundaries in observed energy
  986 + GEbounds ebounds = edisp.ebounds_obs(log10_e_src);
1008 987  
1009 988 // Skip this energy if the boundaries are empty
1010 989 if (ebounds.is_empty()) {
... ... @@ -1012,30 +991,30 @@ void TestGCTAResponse::test_edisp_integration(const GCTAEdisp&amp; edisp,
1012 991 }
1013 992  
1014 993 // Set measured energy range
1015   - GEnergy emin = ebounds.emin();
1016   - GEnergy emax = ebounds.emax();
1017   - double logE_min = std::log10(emin.TeV());
1018   - double logE_max = std::log10(emax.TeV());
1019   -
1020   - // Compute step size for numerical integration
1021   - const int steps = 1000;
1022   - double dlogE = (logE_max-logE_min)/steps;
1023   -
1024   - // Perform numerical integration by summing
1025   - double sum = 0.0;
1026   - double logE_obs = logE_min;
1027   - for (int i = 0; i < steps; ++i) {
1028   - double dp_dlogE = edisp(logE_obs, log10_e_src);
1029   - sum += dp_dlogE * dlogE;
1030   - logE_obs += dlogE;
1031   - }
  994 + GEnergy emin = ebounds.emin();
  995 + GEnergy emax = ebounds.emax();
  996 + double logE_min = std::log10(emin.TeV());
  997 + double logE_max = std::log10(emax.TeV());
  998 +
  999 + // Compute step size for numerical integration
  1000 + const int steps = 1000;
  1001 + double dlogE = (logE_max-logE_min)/steps;
  1002 +
  1003 + // Perform numerical integration by summing
  1004 + double sum = 0.0;
  1005 + double logE_obs = logE_min;
  1006 + for (int i = 0; i < steps; ++i) {
  1007 + double dp_dlogE = edisp(logE_obs, log10_e_src);
  1008 + sum += dp_dlogE * dlogE;
  1009 + logE_obs += dlogE;
  1010 + }
1032 1011  
1033 1012 // Set message string
1034   - GEnergy eng(e_src, "TeV");
  1013 + GEnergy eng(e_src, "TeV");
1035 1014 std::string msg = edisp.classname()+" integration for "+eng.print();
1036   - test_value(sum, 1.0, 0.002, msg);
  1015 + test_value(sum, 1.0, 0.01, msg);
1037 1016  
1038   - }
  1017 + }
1039 1018  
1040 1019 // Return
1041 1020 return;
... ...
src/xspec/GRmf.cpp
  View file @cee203a
... ... @@ -449,9 +449,12 @@ void GRmf::read(const GFitsTable&amp; table)
449 449 // Initialize matrix
450 450 m_matrix = GMatrixSparse(rows, columns);
451 451  
  452 + // Initialize matrix maximum
  453 + double max = 0.0;
  454 +
452 455 // Set true energy bins
453 456 for (int itrue = 0; itrue < rows; ++itrue) {
454   -
  457 +
455 458 // Append energy bin
456 459 GEnergy emin(energy_lo->real(itrue), energy_lo->unit());
457 460 GEnergy emax(energy_hi->real(itrue), energy_hi->unit());
... ... @@ -465,16 +468,24 @@ void GRmf::read(const GFitsTable&amp; table)
465 468 // Get start column index and number of columns
466 469 int imeasured = f_chan->integer(itrue, igroup);
467 470 int nvalues = n_chan->integer(itrue, igroup);
468   -
  471 +
469 472 // Get values
470 473 for (int i = 0; i < nvalues; ++i, ++imeasured, ++icolumn) {
471 474 m_matrix(itrue, imeasured) = matrix->real(itrue, icolumn);
  475 +
  476 + // Update fmax
  477 + if (max < m_matrix(itrue, imeasured)) {
  478 + max = m_matrix(itrue, imeasured);
  479 + m_itruemax = itrue;
  480 + m_imeasmax = imeasured;
  481 + }
472 482 }
473 483  
474 484 } // endfor: looped over groups
475 485  
476 486 } // endfor: looped over true energy bins
477 487  
  488 +
478 489 // Return
479 490 return;
480 491 }
... ... @@ -497,7 +508,7 @@ void GRmf::write(GFits&amp; fits, const std::string&amp; unit) const
497 508 // Compute maximum number of matrix elements
498 509 int maxcolumns = 1;
499 510 for (int i = 0; i < length; i++) {
500   -
  511 +
501 512 // Search first non-zero element
502 513 int ifirst = 0;
503 514 for (; ifirst < nmeasured(); ++ifirst) {
... ... @@ -521,7 +532,7 @@ void GRmf::write(GFits&amp; fits, const std::string&amp; unit) const
521 532 if (num > maxcolumns) {
522 533 maxcolumns = num;
523 534 }
524   -
  535 +
525 536 } // endfor: looped over true energies
526 537  
527 538 // Create new binary table
... ... @@ -560,7 +571,7 @@ void GRmf::write(GFits&amp; fits, const std::string&amp; unit) const
560 571 ifirst = 0;
561 572 num = 0;
562 573 }
563   -
  574 +
564 575 // Set energy bin
565 576 energy_lo(itrue) = m_ebds_true.emin(itrue)(unit);
566 577 energy_hi(itrue) = m_ebds_true.emax(itrue)(unit);
... ... @@ -662,6 +673,8 @@ void GRmf::init_members(void)
662 673 m_ebds_true.clear();
663 674 m_ebds_measured.clear();
664 675 m_matrix.clear();
  676 + m_itruemax = 0;
  677 + m_imeasmax = 0;
665 678  
666 679 // Return
667 680 return;
... ... @@ -680,6 +693,8 @@ void GRmf::copy_members(const GRmf&amp; rmf)
680 693 m_ebds_true = rmf.m_ebds_true;
681 694 m_ebds_measured = rmf.m_ebds_measured;
682 695 m_matrix = rmf.m_matrix;
  696 + m_itruemax = rmf.m_itruemax;
  697 + m_imeasmax = rmf.m_imeasmax;
683 698  
684 699 // Return
685 700 return;
... ...