Justine Devin / gammalib

Commit 5fd8c80d3c9a78c9a0a365d0e95fd81c28fa3881

Authored by Jürgen Knödlseder
1 parent 0c287215

Correct mc() method for the special case that the second and/or third Gaussian are zero.

Inline Side-by-side
Showing 3 changed files with 25 additions and 5 deletions
ChangeLog
  View file @5fd8c80
1   -2014-01-31 Juergen Knoedlseder <jurgen.knodlseder@irap.omp.eu>
  1 +2014-02-04 Juergen Knoedlseder <jurgen.knodlseder@irap.omp.eu>
2 2  
3 3 * Bug fix Version 0.8.1 released
4 4 ================================
... ... @@ -11,6 +11,7 @@
11 11 Added vo module and GVOClient class to Python interface
12 12 Track results of Python unit tests (#1119)
13 13 Clean-up calibration database interface (#1102)
  14 + Correct 2D PSF Monte-Carlo simulations (#1128)
14 15  
15 16  
16 17 2014-01-18 Juergen Knoedlseder <jurgen.knodlseder@irap.omp.eu>
... ...
... ... @@ -2,7 +2,7 @@ New Features and Important Changes in GammaLib 00-08-01
2 2  
3 3 Author(s) : Juergen Knoedlseder
4 4  
5   -31 January 2014
  5 +4 February 2014
6 6  
7 7  
8 8 1. Introduction
... ... @@ -130,6 +130,10 @@ the id attribute to the XML file in case that identifiers are specified.
130 130 Furthermore, GCTAModelRadialAcceptance no longer write the temporal
131 131 model component.
132 132  
  133 +The GCTAPsf2D::mc() method did not sample the PSF correctly in case that
  134 +the sigma and normalization parameters were zero for the second and/or
  135 +third component. This has now been corrected for.
  136 +
133 137  
134 138 21. COMPTEL interface
135 139 ---------------------
... ...
inst/cta/src/GCTAPsf2D.cpp
  View file @5fd8c80
... ... @@ -317,10 +317,10 @@ double GCTAPsf2D::mc(GRan&amp; ran,
317 317 double sum3 = m_sigma3 * m_norm3;
318 318 double sum = sum1 + sum2 + sum3;
319 319 double u = ran.uniform() * sum;
320   - if (u >= sum2) {
  320 + if (sum2 > 0.0 && u >= sum2) {
321 321 sigma = m_sigma3;
322 322 }
323   - else if (u >= sum1) {
  323 + else if (sum1 > 0.0 && u >= sum1) {
324 324 sigma = m_sigma2;
325 325 }
326 326  
... ... @@ -383,12 +383,27 @@ std::string GCTAPsf2D::print(const GChatter&amp; chatter) const
383 383 // Continue only if chatter is not silent
384 384 if (chatter != SILENT) {
385 385  
  386 + // Compute energy boundaries in TeV
  387 + double emin = m_psf.axis_lo(0,0);
  388 + double emax = m_psf.axis_hi(0,m_psf.axis(0)-1);
  389 +
  390 + // Compute offset angle boundaries in deg
  391 + double omin = m_psf.axis_lo(1,0);
  392 + double omax = m_psf.axis_hi(1,m_psf.axis(1)-1);
  393 +
386 394 // Append header
387 395 result.append("=== GCTAPsf2D ===");
388 396  
389 397 // Append information
390 398 result.append("\n"+gammalib::parformat("Filename")+m_filename);
391   - result.append("\n"+m_psf.print(chatter));
  399 + result.append("\n"+gammalib::parformat("Number of energy bins") +
  400 + gammalib::str(m_psf.axis(0)));
  401 + result.append("\n"+gammalib::parformat("Number of offset bins") +
  402 + gammalib::str(m_psf.axis(1)));
  403 + result.append("\n"+gammalib::parformat("Log10(Energy) range"));
  404 + result.append(gammalib::str(emin)+" - "+gammalib::str(emax)+" TeV");
  405 + result.append("\n"+gammalib::parformat("Offset angle range"));
  406 + result.append(gammalib::str(omin)+" - "+gammalib::str(omax)+" deg");
392 407  
393 408 } // endif: chatter was not silent
394 409  
... ...