Olga Sergijenko / gammalib

Commit d6760e1542f2c77f7a66a42506e756af10f9a2bc

Authored by Jürgen Knödlseder
1 parent 5ac42b67

Merge in the current integration branch

Inline Side-by-side
Showing 2 changed files with 69 additions and 43 deletions
include/GResponse.hpp
  View file @d6760e1
... ... @@ -122,14 +122,17 @@ protected:
122 122 void copy_members(const GResponse& rsp);
123 123 void free_members(void);
124 124 double eval_prob(const GModelSky& model,
125   - const GEvent& event,
126   - const GEnergy& srcEng,
127   - const GTime& srcTime,
128   - const GObservation& obs,
129   - const bool& grad) const;
  125 + const GEvent& event,
  126 + const GEnergy& srcEng,
  127 + const GTime& srcTime,
  128 + const GObservation& obs,
  129 + const bool& grad) const;
130 130 GVector eval_probs(const GModelSky& model,
131 131 const GObservation& obs,
132 132 GMatrixSparse* gradients = NULL) const;
  133 + int size_edisp_vector(const GModelSky& model,
  134 + const GObservation& obs,
  135 + const bool& grad) const;
133 136  
134 137 // Virtual protected methods
135 138 virtual double irf_ptsrc(const GEvent& event,
... ...
src/obs/GResponse.cpp
  View file @d6760e1
... ... @@ -207,37 +207,9 @@ double GResponse::convolve(const GModelSky& model,
207 207 GEbounds ebounds = this->ebounds(event.energy());
208 208  
209 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   - }
  210 + int size = size_edisp_vector(model, obs, grad);
239 211  
240   - // Initialise Ndarray array
  212 + // Initialise vector array
241 213 GVector array(size);
242 214  
243 215 // Loop over all boundaries
... ... @@ -410,6 +382,9 @@ GVector GResponse::convolve(const GModelSky&amp; model,
410 382 }
411 383 }
412 384  
  385 + // Determine size of result array
  386 + int size = size_edisp_vector(model, obs, grad);
  387 +
413 388 // Loop over events
414 389 for (int k = 0; k < nevents; ++k) {
415 390  
... ... @@ -422,8 +397,8 @@ GVector GResponse::convolve(const GModelSky&amp; model,
422 397 // Retrieve true energy boundaries
423 398 GEbounds ebounds = this->ebounds(event.energy());
424 399  
425   - // Initialise vector
426   - GVector array;
  400 + // Initialise vector array
  401 + GVector array(size);
427 402  
428 403 // Loop over all boundaries
429 404 for (int i = 0; i < ebounds.size(); ++i) {
... ... @@ -443,12 +418,7 @@ GVector GResponse::convolve(const GModelSky&amp; model,
443 418 integral.fixed_iter(iter);
444 419  
445 420 // Do Romberg integration
446   - if (array.size() == 0) {
447   - array = integral.romberg(etrue_min, etrue_max);
448   - }
449   - else {
450   - array += integral.romberg(etrue_min, etrue_max);
451   - }
  421 + array += integral.romberg(etrue_min, etrue_max);
452 422  
453 423 } // endif: interval was valid
454 424  
... ... @@ -1793,6 +1763,59 @@ GVector GResponse::eval_probs(const GModelSky&amp; model,
1793 1763  
1794 1764  
1795 1765 /***********************************************************************//**
  1766 + * @brief Return size of vector for energy dispersion computation
  1767 + *
  1768 + * @param[in] model Sky model.
  1769 + * @param[in] obs Observation.
  1770 + * @param[out] grad Signals whether gradients should be computed.
  1771 + * @return Size of vector for energy dispersion computation.
  1772 + *
  1773 + * Computes the size of the vector that will be computed for the computation
  1774 + * of the energy dispersion.
  1775 + ***************************************************************************/
  1776 +int GResponse::size_edisp_vector(const GModelSky& model,
  1777 + const GObservation& obs,
  1778 + const bool& grad) const
  1779 +{
  1780 + // Initialise vector size
  1781 + int size = 1;
  1782 +
  1783 + // Determine size of gradient part
  1784 + if (grad) {
  1785 + if (model.spectral() != NULL) {
  1786 + for (int i = 0; i < model.spectral()->size(); ++i) {
  1787 + GModelPar& par = (*(model.spectral()))[i];
  1788 + if (par.is_free() && par.has_grad()) {
  1789 + size++;
  1790 + }
  1791 + }
  1792 + }
  1793 + if (model.temporal() != NULL) {
  1794 + for (int i = 0; i < model.temporal()->size(); ++i) {
  1795 + GModelPar& par = (*(model.temporal()))[i];
  1796 + if (par.is_free() && par.has_grad()) {
  1797 + size++;
  1798 + }
  1799 + }
  1800 + }
  1801 + if (model.has_scales()) {
  1802 + for (int i = 0; i < model.scales(); ++i) {
  1803 + GModelPar& par = const_cast<GModelPar&>(model.scale(i));
  1804 + if (par.name() == obs.instrument()) {
  1805 + if (par.is_free() && par.has_grad()) {
  1806 + size++;
  1807 + }
  1808 + }
  1809 + }
  1810 + }
  1811 + }
  1812 +
  1813 + // Return size
  1814 + return size;
  1815 +}
  1816 +
  1817 +
  1818 +/***********************************************************************//**
1796 1819 * @brief Constructor for energy dispersion integration kernels class
1797 1820 *
1798 1821 * @param[in] parent Pointer to response.
... ...