Commit 504e92677a0febbdb8fae5754ff3e1b40e5c34b4

Authored by Jürgen Knödlseder
1 parent 9bf0d2f7

Improve CTA King PSF code

Showing 1 changed file with 22 additions and 30 deletions
inst/cta/src/GCTAPsfKing.cpp
@@ -201,36 +201,28 @@ double GCTAPsfKing::operator()(const double& delta, @@ -201,36 +201,28 @@ double GCTAPsfKing::operator()(const double& delta,
201 // Update the parameter cache 201 // Update the parameter cache
202 update(logE, theta); 202 update(logE, theta);
203 203
204 - #if defined(G_SMOOTH_PSF)  
205 - // Set ramp down radius  
206 - double ramp_down = 0.95 * m_par_rmax;  
207 - double norm_down = 1.0 / (m_par_rmax - ramp_down);  
208 - #endif  
209 -  
210 - // Continue only if delta is smaller than PSF radius  
211 - if (delta <= m_par_rmax) {  
212 -  
213 - // Continue only if normalization is positive  
214 - if (m_par_norm > 0.0) {  
215 -  
216 - // Compute PSF value  
217 - double arg = delta / m_par_sigma;  
218 - double arg2 = arg * arg;  
219 - psf = m_par_norm *  
220 - std::pow((1.0 + 1.0 / (2.0 * m_par_gamma) * arg2), -m_par_gamma);  
221 -  
222 - // If we are at large offset angles, add a smooth ramp down to  
223 - // avoid steps in the log-likelihood computation  
224 - #if defined(G_SMOOTH_PSF)  
225 - if (delta > ramp_down) {  
226 - double x = norm_down * (delta - ramp_down);  
227 - psf *= 1.0 - x * x;  
228 - }  
229 - #endif  
230 -  
231 - } // endif: normalization was positive  
232 -  
233 - } // endif: PSF radius was smaller than PSF radius 204 + // Continue only if delta is smaller than PSF radius and normalization is
  205 + // positive
  206 + if ((delta <= m_par_rmax) && (m_par_norm > 0.0)) {
  207 +
  208 + // Compute PSF value
  209 + double arg = delta / m_par_sigma;
  210 + double arg2 = arg * arg;
  211 + psf = m_par_norm *
  212 + std::pow((1.0 + 1.0 / (2.0 * m_par_gamma) * arg2), -m_par_gamma);
  213 +
  214 + // If we are at large offset angles, add a smooth ramp down to
  215 + // avoid steps in the log-likelihood computation
  216 + #if defined(G_SMOOTH_PSF)
  217 + double ramp_down = 0.95 * m_par_rmax;
  218 + double norm_down = 1.0 / (m_par_rmax - ramp_down);
  219 + if (delta > ramp_down) {
  220 + double x = norm_down * (delta - ramp_down);
  221 + psf *= 1.0 - x * x;
  222 + }
  223 + #endif
  224 +
  225 + } // endif: PSF radius was smaller than PSF radius and normalization was > 0
234 226
235 // Return PSF 227 // Return PSF
236 return psf; 228 return psf;