Eirik Hatlen / gammalib

Commit 5ac42b678bf81f4fdd6ec72fde5514a095367d34

Authored by Jürgen Knödlseder
2 parents d064542f ee704e03

Merge branch 'integration' into 3203-parallel-spatial-irf-integration 

# Conflicts:
#	src/obs/GObservation.cpp
#	src/obs/GResponse.cpp
Inline Side-by-side
Showing 2 changed files with 47 additions and 21 deletions
src/obs/GObservation.cpp
  View file @5ac42b6
@@ -1184,7 +1184,7 @@ bool GObservation::has_gradient(const GModelPar& par) const @@ -1184,7 +1184,7 @@ bool GObservation::has_gradient(const GModelPar& par) const
1184 { 1184 {
1185 // Initialise flag 1185 // Initialise flag
1186 bool has_gradient = false; 1186 bool has_gradient = false;
1187 - 1187 +
1188 // Only consider free parameters 1188 // Only consider free parameters
1189 if (par.is_free()) { 1189 if (par.is_free()) {
1190 1190
@@ -1374,10 +1374,10 @@ double GObservation::likelihood_poisson_unbinned(const GModels& models, @@ -1374,10 +1374,10 @@ double GObservation::likelihood_poisson_unbinned(const GModels& models,
1374 // Create index array of non-zero derivatives and initialise working 1374 // Create index array of non-zero derivatives and initialise working
1375 // array 1375 // array
1376 int ndev = 0; 1376 int ndev = 0;
1377 - for (int i = 0; i < npars; ++i) {  
1378 - values[i] = 0.0;  
1379 - if (wrk_grad[i] != 0.0 && !gammalib::is_infinite(wrk_grad[i])) {  
1380 - inx[ndev] = i; 1377 + for (int ipar = 0; ipar < npars; ++ipar) {
  1378 + values[ipar] = 0.0;
  1379 + if (wrk_grad[ipar] != 0.0 && !gammalib::is_infinite(wrk_grad[ipar])) {
  1380 + inx[ndev] = ipar;
1381 ndev++; 1381 ndev++;
1382 } 1382 }
1383 } 1383 }
@@ -1538,10 +1538,10 @@ double GObservation::likelihood_poisson_binned(const GModels&amp; models, @@ -1538,10 +1538,10 @@ double GObservation::likelihood_poisson_binned(const GModels&amp; models,
1538 // Create index array of non-zero derivatives and initialise working 1538 // Create index array of non-zero derivatives and initialise working
1539 // array 1539 // array
1540 int ndev = 0; 1540 int ndev = 0;
1541 - for (int k = 0; k < npars; ++k) {  
1542 - values[k] = 0.0;  
1543 - if (wrk_grad[k] != 0.0 && !gammalib::is_infinite(wrk_grad[k])) {  
1544 - inx[ndev] = k; 1541 + for (int ipar = 0; ipar < npars; ++ipar) {
  1542 + values[ipar] = 0.0;
  1543 + if (wrk_grad[ipar] != 0.0 && !gammalib::is_infinite(wrk_grad[ipar])) {
  1544 + inx[ndev] = ipar;
1545 ndev++; 1545 ndev++;
1546 } 1546 }
1547 } 1547 }
@@ -1726,10 +1726,10 @@ double GObservation::likelihood_gaussian_binned(const GModels&amp; models, @@ -1726,10 +1726,10 @@ double GObservation::likelihood_gaussian_binned(const GModels&amp; models,
1726 // Create index array of non-zero derivatives and initialise working 1726 // Create index array of non-zero derivatives and initialise working
1727 // array 1727 // array
1728 int ndev = 0; 1728 int ndev = 0;
1729 - for (int i = 0; i < npars; ++i) {  
1730 - values[i] = 0.0;  
1731 - if (wrk_grad[i] != 0.0 && !gammalib::is_infinite(wrk_grad[i])) {  
1732 - inx[ndev] = i; 1729 + for (int ipar = 0; ipar < npars; ++ipar) {
  1730 + values[ipar] = 0.0;
  1731 + if (wrk_grad[ipar] != 0.0 && !gammalib::is_infinite(wrk_grad[ipar])) {
  1732 + inx[ndev] = ipar;
1733 ndev++; 1733 ndev++;
1734 } 1734 }
1735 } 1735 }
src/obs/GResponse.cpp
  View file @5ac42b6
@@ -206,8 +206,39 @@ double GResponse::convolve(const GModelSky&amp; model, @@ -206,8 +206,39 @@ double GResponse::convolve(const GModelSky&amp; model,
206 // Retrieve true energy boundaries 206 // Retrieve true energy boundaries
207 GEbounds ebounds = this->ebounds(event.energy()); 207 GEbounds ebounds = this->ebounds(event.energy());
208 208
209 - // Initialise vector  
210 - GVector array; 209 + // Determine size of result array
  210 + int size = 1;
  211 + if (grad) {
  212 + if (model.spectral() != NULL) {
  213 + for (int i = 0; i < model.spectral()->size(); ++i) {
  214 + GModelPar& par = (*(model.spectral()))[i];
  215 + if (par.is_free() && par.has_grad()) {
  216 + size++;
  217 + }
  218 + }
  219 + }
  220 + if (model.temporal() != NULL) {
  221 + for (int i = 0; i < model.temporal()->size(); ++i) {
  222 + GModelPar& par = (*(model.temporal()))[i];
  223 + if (par.is_free() && par.has_grad()) {
  224 + size++;
  225 + }
  226 + }
  227 + }
  228 + if (model.has_scales()) {
  229 + for (int i = 0; i < model.scales(); ++i) {
  230 + GModelPar& par = const_cast<GModelPar&>(model.scale(i));
  231 + if (par.name() == obs.instrument()) {
  232 + if (par.is_free() && par.has_grad()) {
  233 + size++;
  234 + }
  235 + }
  236 + }
  237 + }
  238 + }
  239 +
  240 + // Initialise Ndarray array
  241 + GVector array(size);
211 242
212 // Loop over all boundaries 243 // Loop over all boundaries
213 for (int i = 0; i < ebounds.size(); ++i) { 244 for (int i = 0; i < ebounds.size(); ++i) {
@@ -227,12 +258,7 @@ double GResponse::convolve(const GModelSky&amp; model, @@ -227,12 +258,7 @@ double GResponse::convolve(const GModelSky&amp; model,
227 integral.fixed_iter(iter); 258 integral.fixed_iter(iter);
228 259
229 // Do Romberg integration 260 // Do Romberg integration
230 - if (array.size() == 0) {  
231 - array = integral.romberg(etrue_min, etrue_max);  
232 - }  
233 - else {  
234 - array += integral.romberg(etrue_min, etrue_max);  
235 - } 261 + array += integral.romberg(etrue_min, etrue_max);
236 262
237 } // endif: interval was valid 263 } // endif: interval was valid
238 264