Fix energy dispersion bug after implementation of analytical gradients (#3403)

Jürgen Knödlseder
1 parent e0888fb7
Showing 5 changed files with 54 additions and 22 deletions
src/model/GModelSpatialDiffuse.cpp
src/model/GModelSpatialElliptical.cpp
src/model/GModelSpatialRadial.cpp
src/obs/GObservation.cpp
src/obs/GResponse.cpp
@@ -65,7 +65,7 @@ GModelSpatialDiffuse::GModelSpatialDiffuse(void) : GModelSpatial()
  * @param[in] model Diffuse spatial model.
  ***************************************************************************/
 GModelSpatialDiffuse::GModelSpatialDiffuse(const GModelSpatialDiffuse& model) :
-                      GModelSpatial()
+                      GModelSpatial(model)
 {
     // Initialise members
     init_members();
@@ -89,7 +89,7 @@ GModelSpatialElliptical::GModelSpatialElliptical(const GXmlElement&amp; xml) :
  * @param[in] model Elliptical spatial model.
  ***************************************************************************/
 GModelSpatialElliptical::GModelSpatialElliptical(const GModelSpatialElliptical& model) :
-                         GModelSpatial()
+                         GModelSpatial(model)
 {
     // Initialise members
     init_members();
@@ -90,7 +90,7 @@ GModelSpatialRadial::GModelSpatialRadial(const GXmlElement&amp; xml) :
  * @param[in] model Radial spatial model.
  ***************************************************************************/
 GModelSpatialRadial::GModelSpatialRadial(const GModelSpatialRadial& model) :
-                     GModelSpatial()
+                     GModelSpatial(model)
 {
     // Initialise members
     init_members();
@@ -1006,10 +1006,10 @@ double GObservation::likelihood_poisson_unbinned(const GModels&amp; models,
         // Create index array of non-zero derivatives and initialise working
         // array
         int ndev = 0;
-        for (int i = 0; i < npars; ++i) {
-            values[i] = 0.0;
-            if (wrk_grad[i] != 0.0 && !gammalib::is_infinite(wrk_grad[i])) {
-                inx[ndev] = i;
+        for (int ipar = 0; ipar < npars; ++ipar) {
+            values[ipar] = 0.0;
+            if (wrk_grad[ipar] != 0.0 && !gammalib::is_infinite(wrk_grad[ipar])) {
+                inx[ndev] = ipar;
                 ndev++;
             }
         }
@@ -1166,13 +1166,19 @@ double GObservation::likelihood_poisson_binned(const GModels&amp; models,
         // Create index array of non-zero derivatives and initialise working
         // array
         int ndev = 0;
-        for (int i = 0; i < npars; ++i) {
-            values[i] = 0.0;
-            if (wrk_grad[i] != 0.0 && !gammalib::is_infinite(wrk_grad[i])) {
-                inx[ndev] = i;
+        for (int ipar = 0; ipar < npars; ++ipar) {
+            values[ipar] = 0.0;
+            if (wrk_grad[ipar] != 0.0 && !gammalib::is_infinite(wrk_grad[ipar])) {
+                inx[ndev] = ipar;
                 ndev++;
             }
+            /*
+            else if gammalib::is_infinite(wrk_grad[ipar]) {
+                std::cout << "Infinite gradient for ipar=" << ipar << std::endl;
+            }
+            */
         }
+//std::cout << ndev << std::endl;
         // Update gradient vector and curvature matrix. To avoid
         // unneccessary computations we distinguish the case where
@@ -1354,10 +1360,10 @@ double GObservation::likelihood_gaussian_binned(const GModels&amp; models,
         // Create index array of non-zero derivatives and initialise working
         // array
         int ndev = 0;
-        for (int i = 0; i < npars; ++i) {
-            values[i] = 0.0;
-            if (wrk_grad[i] != 0.0 && !gammalib::is_infinite(wrk_grad[i])) {
-                inx[ndev] = i;
+        for (int ipar = 0; ipar < npars; ++ipar) {
+            values[ipar] = 0.0;
+            if (wrk_grad[ipar] != 0.0 && !gammalib::is_infinite(wrk_grad[ipar])) {
+                inx[ndev] = ipar;
                 ndev++;
             }
         }
@@ -203,8 +203,39 @@ double GResponse::convolve(const GModelSky&amp;    model,
             // Retrieve true energy boundaries
             GEbounds ebounds = this->ebounds(event.energy());
+            // Determine size of result array
+            int size = 1;
+            if (grad) {
+                if (model.spectral() != NULL) {
+                    for (int i = 0; i < model.spectral()->size(); ++i) {
+                        GModelPar& par = (*(model.spectral()))[i];
+                        if (par.is_free() && par.has_grad()) {
+                            size++;
+                        }
+                    }
+                }
+                if (model.temporal() != NULL) {
+                    for (int i = 0; i < model.temporal()->size(); ++i) {
+                        GModelPar& par = (*(model.temporal()))[i];
+                        if (par.is_free() && par.has_grad()) {
+                            size++;
+                        }
+                    }
+                }
+                if (model.has_scales()) {
+                    for (int i = 0; i < model.scales(); ++i) {
+                        GModelPar& par = const_cast<GModelPar&>(model.scale(i));
+                        if (par.name() == obs.instrument()) {
+                            if (par.is_free() && par.has_grad()) {
+                                size++;
+                            }
+                        }
+                    }
+                }
+            }
+
             // Initialise Ndarray array
-            GNdarray array;
+            GNdarray array(size);
             // Loop over all boundaries
             for (int i = 0; i < ebounds.size(); ++i) {
@@ -224,12 +255,7 @@ double GResponse::convolve(const GModelSky&amp;    model,
                     integral.fixed_iter(iter);
                     // Do Romberg integration
-                    if (array.size() == 0) {
-                        array = integral.romberg(etrue_min, etrue_max);
-                    }
-                    else {
-                        array += integral.romberg(etrue_min, etrue_max);
-                    }
+                    array += integral.romberg(etrue_min, etrue_max);
                 } // endif: interval was valid