Integrate code

Jürgen Knödlseder
1 parent cee203af
Showing 5 changed files with 139 additions and 143 deletions
include/GRmf.hpp
inst/cta/include/GCTAEdisp2D.hpp
inst/cta/include/GCTAEdispRmf.hpp
inst/cta/src/GCTAEdisp2D.cpp
inst/cta/src/GCTAEdispRmf.cpp
@@ -65,6 +65,8 @@ public:
     int                  size(void) const;
     int                  ntrue(void) const;
     int                  nmeasured(void) const;
+    int                  itruemax(void) const;
+    int                  imeasmax(void) const;
     double&              at(const int& itrue, const int& imeasured);
     const double&        at(const int& itrue, const int& imeasured) const;
     const GEbounds&      etrue(void) const;
@@ -81,8 +83,6 @@ public:
                                const std::string& unit = "keV") const;
     const std::string&   filename(void) const;
     std::string          print(const GChatter& chatter = NORMAL) const;
-    int                  itruemax(void) const;
-    int                  imeasmax(void) const;
  
 protected:
     // Protected methods
@@ -95,8 +95,8 @@ protected:
     GEbounds            m_ebds_true;     //!< True energy boundaries
     GEbounds            m_ebds_measured; //!< Measured energy boundaries
     GMatrixSparse       m_matrix;        //!< Sparse redistribution matrix
-    mutable int         m_imeasmax;
-    mutable int         m_itruemax;
+    int                 m_imeasmax;      //!< Index of measured maximum
+    int                 m_itruemax;      //!< Index of true maximum
  
 };
  
@@ -187,6 +187,36 @@ int GRmf::nmeasured(void) const
 }
  
  
+
+
+/***********************************************************************//**
+ * @brief Return true energy index of maximum value of the redistribution matrix
+ *
+ * @return True energy index of maximum value of the redistribution matrix.
+ *
+ * Returns the true energy index of maximum value of the redistribution matrix.
+ ***************************************************************************/
+inline
+int GRmf::itruemax(void) const
+{
+    return m_itruemax;
+}
+
+
+/***********************************************************************//**
+ * @brief Return measured energy index of maximum value of the redistribution matrix
+ *
+ * @return Measured energy index of maximum value of the redistribution matrix.
+ *
+ * Returns the measured energy index of maximum value of the redistribution matrix.
+ ***************************************************************************/
+inline
+int GRmf::imeasmax(void) const
+{
+    return m_imeasmax;
+}
+
+
 /***********************************************************************//**
  * @brief Return true energy boundaries
  *
@@ -243,31 +273,4 @@ const std::string&amp; GRmf::filename(void) const
     return (m_filename);
 }
  
-
-/***********************************************************************//**
- * @brief Return true energy index of maximum value of the redistribution matrix
- *
- * @return True energy index of maximum value of the redistribution matrix.
- *
- * Returns the true energy index of maximum value of the redistribution matrix.
- ***************************************************************************/
-inline
-int GRmf::itruemax(void) const
-{
-    return m_itruemax;
-}
-
-/***********************************************************************//**
- * @brief Return measured energy index of maximum value of the redistribution matrix
- *
- * @return Measured energy index of maximum value of the redistribution matrix.
- *
- * Returns the measured energy index of maximum value of the redistribution matrix.
- ***************************************************************************/
-inline
-int GRmf::imeasmax(void) const
-{
-    return m_imeasmax;
-}
-
 #endif /* GRMF_HPP */
@@ -121,6 +121,7 @@ private:
     std::string       m_filename;  //!< Name of Edisp response file
     GCTAResponseTable m_edisp;     //!< Edisp response table
  
+    // Monte Carlo cache
     mutable bool      m_cdf_computed;
     mutable double    m_theta;
     mutable bool      m_ebounds_obs_computed;
@@ -131,11 +132,11 @@ private:
     mutable int       m_index_src;
     mutable double    m_fmax;
  
-    mutable std::vector<GEbounds> m_ebounds_obs;
-    mutable std::vector<GEbounds> m_ebounds_src;
+    mutable std::vector<GEbounds>                                 m_ebounds_obs;
+    mutable std::vector<GEbounds>                                 m_ebounds_src;
     mutable std::vector<std::vector<std::pair<double, double> > > m_cumul;
-    mutable std::vector<double> m_eobs_axis;
-    mutable std::vector<std::vector<double> > m_temp;
+    mutable std::vector<double>                                   m_eobs_axis;
+    mutable std::vector<std::vector<double> >                     m_temp;
 };
  
  
@@ -141,12 +141,12 @@ private:
     mutable double                                  m_logEobs;
     mutable std::vector<std::vector<std::pair<double, double> > > m_cumul;
  
-    mutable bool m_ebounds_obs_computed;
+    mutable bool                  m_ebounds_obs_computed;
     mutable std::vector<GEbounds> m_ebounds_obs;
-    mutable int m_index_obs;
-    mutable bool m_ebounds_src_computed;
+    mutable int                   m_index_obs;
+    mutable bool                  m_ebounds_src_computed;
     mutable std::vector<GEbounds> m_ebounds_src;
-    mutable int m_index_src;
+    mutable int                   m_index_src;
  
 };
  
@@ -327,25 +327,19 @@ void GCTAEdisp2D::read(const GFits&amp; fits)
     // Set offset angle axis to radians
     m_edisp.axis_radians(2);
  
-
-//    for (int i_theta = 0; i_theta < theta_size; ++i_theta) {
-
-        m_temp.clear();
-
-        int step = migra_size / etrue_size;
-
-        for (int itrue = 0; itrue < etrue_size; ++itrue) {
-            double logEsrc = std::log10(0.5 * (m_edisp.axis_hi(0, itrue) +
-                                               m_edisp.axis_lo(0, itrue)));
-            m_temp.push_back(std::vector<double>());
-            for (int iobs = 0; iobs < etrue_size; ++iobs) {
-                double EobsOverEsrc = 0.5 * (m_edisp.axis_hi(1, iobs*step) +
-                                             m_edisp.axis_lo(1, iobs*step));
-                m_temp[itrue].push_back(operator()(logEsrc, EobsOverEsrc));
-            }
+    // ...
+    m_temp.clear();
+    int step = migra_size / etrue_size;
+    for (int itrue = 0; itrue < etrue_size; ++itrue) {
+        double logEsrc = std::log10(0.5 * (m_edisp.axis_hi(0, itrue) +
+                                           m_edisp.axis_lo(0, itrue)));
+        m_temp.push_back(std::vector<double>());
+        for (int iobs = 0; iobs < etrue_size; ++iobs) {
+            double EobsOverEsrc = 0.5 * (m_edisp.axis_hi(1, iobs*step) +
+                                         m_edisp.axis_lo(1, iobs*step));
+            m_temp[itrue].push_back(operator()(logEsrc, EobsOverEsrc));
         }
-
- //   }
+    }
  
     // Return
     return;
@@ -395,56 +389,6 @@ void GCTAEdisp2D::save(const std::string&amp; filename, const bool&amp; clobber) const
     return;
 }
  
-/***********************************************************************//**
- * @brief Find maximum energy dispersion
- *
-***************************************************************************/
-void GCTAEdisp2D::set_fmax(void) const
-{
-
-    double max = 0.0;
-
-    // Loop over Esrc
-    for (int i = 0; i < m_edisp.axis(0); ++i) {
-
-        // Set Esrc
-        double Esrc   = 0.5 * (m_edisp.axis_hi(0, i) +
-                               m_edisp.axis_lo(0, i));
-        double logEsrc = std::log10(Esrc);
-
-        // Loop over EobsOverEsrc
-        for (int j = 0; j < m_edisp.axis(1); ++j) {
-
-            // Set EobsOverEsrc
-            double EobsOverEsrc = 0.5 * (m_edisp.axis_hi(1, j) +
-                                         m_edisp.axis_lo(1, j));
-/*
-            // Loop over theta
-            for (int k = 0; k < m_edisp.axis(2); ++k) {
-
-                // Set theta
-                double theta = 0.5 * (m_edisp.axis_hi(2, k) +
-                                     m_edisp.axis_lo(2, k));
-
-                double test = operator()(logEsrc, EobsOverEsrc, theta);
-*/
-                double test = operator()(logEsrc, EobsOverEsrc);
-                if (max < test) {
-                    max = test;
-                }
-/*
-             } // endfor: looped over theta
-*/
-        } // endfor: looped over EobsOverEsrc
-
-    } // endfor: looped over Esrc
-
-    // Set fmax
-    m_fmax = max;
-
-    // Return
-    return;
-}
  
 /***********************************************************************//**
  * @brief Simulate energy dispersion
@@ -479,7 +423,6 @@ GEnergy GCTAEdisp2D::mc(GRan&amp;         ran,
     if (m_fmax == -1.0) {
         set_fmax();
     }
-
     double fmax = m_fmax;
  
     // Find energy by rejection method
@@ -493,9 +436,9 @@ GEnergy GCTAEdisp2D::mc(GRan&amp;         ran,
         ftest  = ran.uniform() * fmax;
     }
  
+    // Set energy
     energy.log10TeV(e);
  
-
     // Return energy
     return energy;
  
@@ -760,23 +703,22 @@ void GCTAEdisp2D::init_members(void)
     m_edisp.clear();
  
     // Initialise Monte Carlo cache
-    m_theta        = 0.0;
-    m_logEsrc      =-30.0;
-    m_logEobs      =-30.0;
-    m_fmax         = -1.0;
-    m_index_obs    = 0;
-    m_index_src    = 0;
-    m_cdf_computed = false;
+    m_theta                =   0.0;
+    m_logEsrc              = -30.0;
+    m_logEobs              = -30.0;
+    m_fmax                 =  -1.0;
+    m_index_obs            =     0;
+    m_index_src            =     0;
+    m_cdf_computed         = false;
     m_ebounds_obs_computed = false;
     m_ebounds_src_computed = false;
  
-    m_cumul.clear();
-
-    // Initialize ebounds
+    // Initialize vectors
     m_ebounds_obs.clear();
     m_ebounds_src.clear();
-
+    m_cumul.clear();
     m_eobs_axis.clear();
+    m_temp.clear();
  
     // Return
     return;
@@ -793,17 +735,22 @@ void GCTAEdisp2D::copy_members(const GCTAEdisp2D&amp; edisp)
     // Copy members
     m_filename  = edisp.m_filename;
     m_edisp     = edisp.m_edisp;
-    m_theta     = edisp.m_theta;
-    m_logEsrc   = edisp.m_logEsrc;
-    m_logEobs   = edisp.m_logEobs;
-    m_index_obs = edisp.m_index_obs;
-    m_index_src = edisp.m_index_src;
-
-    m_ebounds_obs = edisp.m_ebounds_obs;
-    m_ebounds_src = edisp.m_ebounds_src;
-    m_cumul       = edisp.m_cumul;
-    m_eobs_axis   = edisp.m_eobs_axis;
-    m_temp        = edisp.m_temp;
+
+    // Copy Monte Carlo cache
+    m_theta                = edisp.m_theta;
+    m_logEsrc              = edisp.m_logEsrc;
+    m_logEobs              = edisp.m_logEobs;
+    m_fmax                 = edisp.m_fmax;
+    m_index_obs            = edisp.m_index_obs;
+    m_index_src            = edisp.m_index_src;
+    m_cdf_computed         = edisp.m_cdf_computed;
+    m_ebounds_obs_computed = edisp.m_ebounds_obs_computed;
+    m_ebounds_src_computed = edisp.m_ebounds_src_computed;
+    m_ebounds_obs          = edisp.m_ebounds_obs;
+    m_ebounds_src          = edisp.m_ebounds_src;
+    m_cumul                = edisp.m_cumul;
+    m_eobs_axis            = edisp.m_eobs_axis;
+    m_temp                 = edisp.m_temp;
  
     // Return
     return;
@@ -1118,3 +1065,41 @@ void GCTAEdisp2D::compute_cumul(const double&amp; theta,
     return;
 }
 #endif
+
+
+/***********************************************************************//**
+ * @brief Find maximum energy dispersion
+ ***************************************************************************/
+void GCTAEdisp2D::set_fmax(void) const
+{
+    // Initialise maximum
+    double max = 0.0;
+
+    // Loop over Esrc
+    for (int i = 0; i < m_edisp.axis(0); ++i) {
+
+        // Set Esrc
+        double Esrc    = 0.5 * (m_edisp.axis_hi(0,i) + m_edisp.axis_lo(0,i));
+        double logEsrc = std::log10(Esrc);
+
+        // Loop over EobsOverEsrc
+        for (int j = 0; j < m_edisp.axis(1); ++j) {
+
+            // Set EobsOverEsrc
+            double EobsOverEsrc = 0.5 * (m_edisp.axis_hi(1, j) +
+                                         m_edisp.axis_lo(1, j));
+            double test = operator()(logEsrc, EobsOverEsrc);
+            if (max < test) {
+                max = test;
+            }
+
+        } // endfor: looped over EobsOverEsrc
+
+    } // endfor: looped over Esrc
+
+    // Set fmax
+    m_fmax = max;
+
+    // Return
+    return;
+}
@@ -715,18 +715,17 @@ void GCTAEdispRmf::init_members(void)
     // Initialise Monte Carlo cache
     m_mc_measured_start.clear();
     m_mc_measured_cdf.clear();
-    m_cdf_computed   = false;
-    m_theta          = 0.0;
-    m_logEsrc        = -30.0;
-    m_logEobs        = -30.0;
+    m_cdf_computed         = false;
+    m_theta                = 0.0;
+    m_logEsrc              = -30.0;
+    m_logEobs              = -30.0;
     m_cumul.clear();
-
     m_ebounds_obs_computed = false;
     m_ebounds_obs.clear();
-    m_index_obs = 0;
+    m_index_obs            = 0;
     m_ebounds_src_computed = false;
     m_ebounds_src.clear();
-    m_index_src = 0;
+    m_index_src            = 0;
  
     // Return
     return;
@@ -760,11 +759,19 @@ void GCTAEdispRmf::copy_members(const GCTAEdispRmf&amp; edisp)
     m_wgt4           = edisp.m_wgt4;
  
     // Copy Monte Carlo cache
-    m_mc_measured_start = edisp.m_mc_measured_start;
-    m_mc_measured_cdf   = edisp.m_mc_measured_cdf;
-    m_theta             = edisp.m_theta;
-    m_cumul             = edisp.m_cumul;
-    m_cdf_computed      = edisp.m_cdf_computed;
+    m_mc_measured_start    = edisp.m_mc_measured_start;
+    m_mc_measured_cdf      = edisp.m_mc_measured_cdf;
+    m_cdf_computed         = edisp.m_cdf_computed;
+    m_theta                = edisp.m_theta;
+    m_logEsrc              = edisp.m_logEsrc;
+    m_logEobs              = edisp.m_logEobs;
+    m_cumul                = edisp.m_cumul;
+    m_ebounds_obs_computed = edisp.m_ebounds_obs_computed;
+    m_ebounds_obs          = edisp.m_ebounds_obs;
+    m_index_obs            = edisp.m_index_obs;
+    m_ebounds_src_computed = edisp.m_ebounds_src_computed;
+    m_ebounds_src          = edisp.m_ebounds_src;
+    m_index_src            = edisp.m_index_src;
  
     // Return
     return;