Rolf Buehler / gammalib

Commit 2ac11cb013d5c9cc90602de7bb820cc4398f0fe4

Authored by Jürgen Knödlseder
1 parent cc6116bd

Integrate containment_radius code.

Inline Side-by-side
Showing 10 changed files with 315 additions and 282 deletions
inst/cta/include/GCTAPsf.hpp
  View file @2ac11cb
1 1 /***************************************************************************
2 2 * GCTAPsf.hpp - CTA point spread function base class *
3 3 * ----------------------------------------------------------------------- *
4   - * copyright (C) 2012-2013 by Juergen Knoedlseder *
  4 + * copyright (C) 2012-2015 by Juergen Knoedlseder *
5 5 * ----------------------------------------------------------------------- *
6 6 * *
7 7 * This program is free software: you can redistribute it and/or modify *
... ... @@ -81,13 +81,13 @@ public:
81 81 const double& zenith = 0.0,
82 82 const double& azimuth = 0.0,
83 83 const bool& etrue = true) const = 0;
84   - virtual double containment_radius( const double& fraction,
85   - const double& logE,
86   - const double& theta = 0.0,
87   - const double& phi = 0.0,
88   - const double& zenith = 0.0,
89   - const double& azimuth = 0.0,
90   - const bool & etrue = true) const = 0;
  84 + virtual double containment_radius(const double& fraction,
  85 + const double& logE,
  86 + const double& theta = 0.0,
  87 + const double& phi = 0.0,
  88 + const double& zenith = 0.0,
  89 + const double& azimuth = 0.0,
  90 + const bool & etrue = true) const = 0;
91 91 virtual std::string print(const GChatter& chatter = NORMAL) const = 0;
92 92  
93 93 protected:
... ...
inst/cta/include/GCTAPsf2D.hpp
  View file @2ac11cb
1 1 /***************************************************************************
2 2 * GCTAPsf2D.hpp - CTA 2D point spread function class *
3 3 * ----------------------------------------------------------------------- *
4   - * copyright (C) 2012-2014 by Juergen Knoedlseder *
  4 + * copyright (C) 2012-2015 by Juergen Knoedlseder *
5 5 * ----------------------------------------------------------------------- *
6 6 * *
7 7 * This program is free software: you can redistribute it and/or modify *
... ...
inst/cta/include/GCTAPsfKing.hpp
  View file @2ac11cb
1 1 /***************************************************************************
2 2 * GCTAPsfKing.hpp - King profile CTA point spread function class *
3 3 * ----------------------------------------------------------------------- *
4   - * copyright (C) 2013-2014 by Michael Mayer *
  4 + * copyright (C) 2013-2015 by Michael Mayer *
5 5 * ----------------------------------------------------------------------- *
6 6 * *
7 7 * This program is free software: you can redistribute it and/or modify *
... ...
inst/cta/include/GCTAPsfPerfTable.hpp
  View file @2ac11cb
1 1 /***************************************************************************
2 2 * GCTAPsfPerfTable.hpp - CTA performance table PSF class *
3 3 * ----------------------------------------------------------------------- *
4   - * copyright (C) 2012-2014 by Juergen Knoedlseder *
  4 + * copyright (C) 2012-2015 by Juergen Knoedlseder *
5 5 * ----------------------------------------------------------------------- *
6 6 * *
7 7 * This program is free software: you can redistribute it and/or modify *
... ... @@ -83,13 +83,13 @@ public:
83 83 const double& zenith = 0.0,
84 84 const double& azimuth = 0.0,
85 85 const bool& etrue = true) const;
86   - double containment_radius(const double& fraction,
87   - const double& logE,
88   - const double& theta = 0.0,
89   - const double& phi = 0.0,
90   - const double& zenith = 0.0,
91   - const double& azimuth = 0.0,
92   - const bool& etrue = true) const;
  86 + double containment_radius(const double& fraction,
  87 + const double& logE,
  88 + const double& theta = 0.0,
  89 + const double& phi = 0.0,
  90 + const double& zenith = 0.0,
  91 + const double& azimuth = 0.0,
  92 + const bool& etrue = true) const;
93 93 std::string print(const GChatter& chatter = NORMAL) const;
94 94  
95 95 private:
... ...
inst/cta/include/GCTAPsfVector.hpp
  View file @2ac11cb
1 1 /***************************************************************************
2 2 * GCTAPsfVector.hpp - CTA point spread function vector class *
3 3 * ----------------------------------------------------------------------- *
4   - * copyright (C) 2012-2014 by Juergen Knoedlseder *
  4 + * copyright (C) 2012-2015 by Juergen Knoedlseder *
5 5 * ----------------------------------------------------------------------- *
6 6 * *
7 7 * This program is free software: you can redistribute it and/or modify *
... ...
inst/cta/src/GCTAPsf2D.cpp
  View file @2ac11cb
... ... @@ -40,6 +40,8 @@
40 40 /* __ Method name definitions ____________________________________________ */
41 41 #define G_READ "GCTAPsf2D::read(GFits&)"
42 42 #define G_LOAD "GCTAPsf2D::load(std::string&)"
  43 +#define G_CONTAINMENT_RADIUS "GCTAPsf2D::containment_radius(double&,"\
  44 + " double&, double&, double&, double&, double&, bool&)"
43 45  
44 46 /* __ Macros _____________________________________________________________ */
45 47  
... ... @@ -493,6 +495,110 @@ double GCTAPsf2D::delta_max(const double& logE,
493 495  
494 496  
495 497 /***********************************************************************//**
  498 + * @brief Return the radius that contains a fraction of the events (radians)
  499 + *
  500 + * @param[in] fraction of events (0.0<fraction<1.0)
  501 + * @param[in] logE Log10 of the true photon energy (TeV).
  502 + * @param[in] theta Offset angle in camera system (rad). Not used.
  503 + * @param[in] phi Azimuth angle in camera system (rad). Not used.
  504 + * @param[in] zenith Zenith angle in Earth system (rad). Not used.
  505 + * @param[in] azimuth Azimuth angle in Earth system (rad). Not used.
  506 + * @param[in] etrue Use true energy (true/false). Not used.
  507 + * @return Containment radius (radians).
  508 + *
  509 + * @exception GException::invalid_argument
  510 + * Invalid fraction specified.
  511 + *
  512 + * Uses the Newton-Raphson method to find which 'a' solves the equation:
  513 + *
  514 + * \f[
  515 + *
  516 + * fraction = \pi * m\_norm * \left(
  517 + * \frac{ 1 }{m\_width1} e^{m\_width1 * a^2} +
  518 + * \frac{m\_norm2}{m\_width2} e^{m\_width2 * a^2} +
  519 + * \frac{m\_norm3}{m\_width3} e^{m\_width3 * a^2} \right)
  520 + *
  521 + * \f]
  522 + *
  523 + * Calculate the radius from the center that contains 'fraction' percent
  524 + * of the events. fraction * 100. = Containment % . Fraction must be
  525 + * 0.0 < fraction < 1.0 .
  526 + ***************************************************************************/
  527 +double GCTAPsf2D::containment_radius(const double& fraction,
  528 + const double& logE,
  529 + const double& theta,
  530 + const double& phi,
  531 + const double& zenith,
  532 + const double& azimuth,
  533 + const bool& etrue) const
  534 +{
  535 + // Check input argument
  536 + if (fraction <= 0.0 || fraction >= 1.0) {
  537 + std::string message = "Containment fraction "+
  538 + gammalib::str(fraction)+" must be between " +
  539 + "0.0 and 1.0, not inclusive.";
  540 + throw GException::invalid_argument(G_CONTAINMENT_RADIUS, message);
  541 + }
  542 +
  543 + // Update the parameter cache
  544 + update(logE, theta);
  545 +
  546 + // Required accuracy
  547 + const double convergence = 1.0e-6;
  548 +
  549 + // Initial radius to start Newton's method with guess fraction * delta_max
  550 + double a = fraction * delta_max(logE, theta, phi, zenith, azimuth, etrue);
  551 +
  552 + // Maximum number of newton-raphson loops before giving up
  553 + const int iterlimit = 10000;
  554 +
  555 + // Do the newton-raphson loops
  556 + int iter(0);
  557 + for (; iter < iterlimit; ++iter) {
  558 +
  559 + // ...
  560 + double a2 = a*a;
  561 +
  562 + // Calculate f(a)
  563 + double fa(0.0);
  564 + fa += m_norm * (std::exp( m_width1 * a2) - 1.0) / m_width1;
  565 + fa += m_norm2 * (std::exp( m_width2 * a2) - 1.0) / m_width2;
  566 + fa += m_norm3 * (std::exp( m_width3 * a2) - 1.0) / m_width3;
  567 + fa *= gammalib::pi;
  568 + fa -= fraction;
  569 +
  570 + // Check if we've met the desired accuracy
  571 + if (std::abs(fa) < convergence) {
  572 + break;
  573 + }
  574 +
  575 + // Calculate f'(a)
  576 + double fp(0.0);
  577 + fp += m_norm * std::exp(m_width1 * a2);
  578 + fp += m_norm2 * std::exp(m_width2 * a2);
  579 + fp += m_norm3 * std::exp(m_width3 * a2);
  580 + fp *= gammalib::twopi * a;
  581 +
  582 + // Calculate next point via x+1 = x - f(a) / f'(a)
  583 + a = a - fa / fp;
  584 +
  585 + } // endfor: Newton-Raphson loops
  586 +
  587 + // Warn the user if we didn't converge
  588 + if (iter == iterlimit-1) {
  589 + std::string message = "Unable to converge within " +
  590 + gammalib::str(convergence) +
  591 + " of the root in less than " +
  592 + gammalib::str(iterlimit) + " iterations." ;
  593 + gammalib::warning(G_CONTAINMENT_RADIUS, message);
  594 + }
  595 +
  596 + // Return containment radius
  597 + return a;
  598 +}
  599 +
  600 +
  601 +/***********************************************************************//**
496 602 * @brief Print point spread function information
497 603 *
498 604 * @param[in] chatter Chattiness (defaults to NORMAL).
... ... @@ -670,112 +776,3 @@ void GCTAPsf2D::update(const double&amp; logE, const double&amp; theta) const
670 776 // Return
671 777 return;
672 778 }
673   -
674   -
675   -/***********************************************************************//**
676   - * @brief Return the radius that contains a fraction of the events (radians)
677   - *
678   - * @param[in] fraction of events (0.0<fraction<1.0)
679   - * @param[in] logE Log10 of the true photon energy (TeV).
680   - * @param[in] theta Offset angle in camera system (rad). Not used.
681   - * @param[in] phi Azimuth angle in camera system (rad). Not used.
682   - * @param[in] zenith Zenith angle in Earth system (rad). Not used.
683   - * @param[in] azimuth Azimuth angle in Earth system (rad). Not used.
684   - * @param[in] etrue Use true energy (true/false). Not used.
685   - *
686   - * Uses the Newton-Raphson method to find which 'a' solves the equation:
687   - *
688   - * \f[
689   - *
690   - * fraction = \pi * m\_norm * \left(
691   - * \frac{ 1 }{m\_width1} e^{m\_width1 * a^2} +
692   - * \frac{m\_norm2}{m\_width2} e^{m\_width2 * a^2} +
693   - * \frac{m\_norm3}{m\_width3} e^{m\_width3 * a^2} \right)
694   - *
695   - * \f]
696   - *
697   - * Calculate the radius from the center that contains 'fraction' percent
698   - * of the events. fraction * 100. = Containment % . Fraction must be
699   - * 0.0 < fraction < 1.0 .
700   - ***************************************************************************/
701   -double GCTAPsf2D::containment_radius(const double& fraction,
702   - const double& logE,
703   - const double& theta,
704   - const double& phi,
705   - const double& zenith,
706   - const double& azimuth,
707   - const bool& etrue) const
708   -{
709   -
710   - if ( fraction <= 0.0 || fraction >= 1.0 ) {
711   - std::string origin = "GCTAPsf2D::containment_radius" ;
712   - std::string message = " First argument (fraction=" +
713   - gammalib::str(fraction) + ") must be between " +
714   - "0.0 and 1.0, not inclusive." ;
715   - throw GException::invalid_argument( origin, message ) ;
716   - }
717   -
718   -
719   - // Update the parameter cache
720   - update(logE, theta);
721   -
722   - // required accuracy
723   - double convergence = 1.0e-6 ;
724   -
725   - // initial radius to start Newton's method with
726   - // guess fraction * delta_max
727   - double a = fraction * delta_max( logE, theta, phi, zenith, azimuth, etrue ) ;
728   -
729   - // maximum number of newton-raphson loops before giving up
730   - int iterlimit = 10000 ;
731   -
732   - // initialization variables for f(a), f'(a), a^2, and the iter counter
733   - double fa = 0.0 ;
734   - double fp = 0.0 ;
735   - double a2 = 0.0 ;
736   - int iter = 0 ;
737   -
738   - // do the newton-raphson loops
739   - for ( iter = 0 ; iter < iterlimit ; iter++ ) {
740   -
741   - a2 = a*a ;
742   -
743   - // calculate f(a)
744   - fa = 0.0 ;
745   - fa += m_norm * ( std::exp( m_width1 * a2 ) - 1 ) / m_width1 ;
746   - fa += m_norm2 * ( std::exp( m_width2 * a2 ) - 1 ) / m_width2 ;
747   - fa += m_norm3 * ( std::exp( m_width3 * a2 ) - 1 ) / m_width3 ;
748   - fa *= gammalib::pi ;
749   - fa -= fraction ;
750   -
751   - // check if we've met the desired accuracy
752   - if ( fabs( fa ) < convergence ) {
753   - break ;
754   - }
755   -
756   - // calculate f'(a)
757   - fp = 0.0 ;
758   - fp += m_norm * std::exp( m_width1 * a2 ) ;
759   - fp += m_norm2 * std::exp( m_width2 * a2 ) ;
760   - fp += m_norm3 * std::exp( m_width3 * a2 ) ;
761   - fp *= gammalib::twopi * a ;
762   -
763   - // calculate next point via x+1 = x - f(a) / f'(a)
764   - a = a - fa / fp ;
765   -
766   - } ;
767   -
768   - // warn the user if we didn't converge
769   - if ( iter == iterlimit-1 ) {
770   - std::string origin = "GCTAPsf2D::containment_radius" ;
771   - std::string message = " Unable to converge within " +
772   - gammalib::str(convergence) +
773   - " of the root in less than " +
774   - gammalib::str(iterlimit) + " iterations." ;
775   - gammalib::warning(origin, message);
776   - }
777   -
778   - // Return containment radius
779   - return a ;
780   -}
781   -
... ...
inst/cta/src/GCTAPsfKing.cpp
  View file @2ac11cb
... ... @@ -39,6 +39,8 @@
39 39  
40 40 /* __ Method name definitions ____________________________________________ */
41 41 #define G_READ "GCTAPsfKing::read(GFits&)"
  42 +#define G_CONTAINMENT_RADIUS "GCTAPsfKing::containment_radius(double&,"\
  43 + " double&, double&, double&, double&, double&, bool&)"
42 44 #define G_UPDATE "GCTAPsfKing::update(double&, double&)"
43 45  
44 46 /* __ Macros _____________________________________________________________ */
... ... @@ -493,6 +495,52 @@ double GCTAPsfKing::delta_max(const double&amp; logE,
493 495  
494 496  
495 497 /***********************************************************************//**
  498 + * @brief Return the radius that contains a fraction of the events (radians)
  499 + *
  500 + * @param[in] fraction of events (0.0-1.0)
  501 + * @param[in] logE Log10 of the true photon energy (TeV).
  502 + * @param[in] theta Offset angle in camera system (rad). Not used.
  503 + * @param[in] phi Azimuth angle in camera system (rad). Not used.
  504 + * @param[in] zenith Zenith angle in Earth system (rad). Not used.
  505 + * @param[in] azimuth Azimuth angle in Earth system (rad). Not used.
  506 + * @param[in] etrue Use true energy (true/false). Not used.
  507 + * @return Containment radius (radians).
  508 + *
  509 + * @exception GException::invalid_argument
  510 + * Invalid fraction specified.
  511 + *
  512 + * Calculate the radius from the center that contains 'fraction' percent
  513 + * of the events. fraction * 100. = Containment %.
  514 + ***************************************************************************/
  515 +double GCTAPsfKing::containment_radius(const double& fraction,
  516 + const double& logE,
  517 + const double& theta,
  518 + const double& phi,
  519 + const double& zenith,
  520 + const double& azimuth,
  521 + const bool& etrue) const
  522 +{
  523 + // Check input argument
  524 + if (fraction <= 0.0 || fraction >= 1.0) {
  525 + std::string message = "Containment fraction "+
  526 + gammalib::str(fraction)+" must be between " +
  527 + "0.0 and 1.0, not inclusive.";
  528 + throw GException::invalid_argument(G_CONTAINMENT_RADIUS, message);
  529 + }
  530 +
  531 + // Update the parameter cache
  532 + update(logE, theta);
  533 +
  534 + // Use analytic calculation
  535 + double arg = std::pow(1.0 - fraction, 1.0/(1.0-m_par_gamma));
  536 + double radius = m_par_sigma * std::pow(2.0 * m_par_gamma * (arg - 1.0), 0.5);
  537 +
  538 + // Return radius containing fraction of events
  539 + return radius;
  540 +}
  541 +
  542 +
  543 +/***********************************************************************//**
496 544 * @brief Print point spread function information
497 545 *
498 546 * @return Content of point spread function instance.
... ... @@ -644,42 +692,3 @@ void GCTAPsfKing::update(const double&amp; logE, const double&amp; theta) const
644 692 // Return
645 693 return;
646 694 }
647   -
648   -
649   -/***********************************************************************//**
650   - * @brief Return the radius that contains a fraction of the events (radians)
651   - *
652   - * @param[in] fraction of events (0.0-1.0)
653   - * @param[in] logE Log10 of the true photon energy (TeV).
654   - * @param[in] theta Offset angle in camera system (rad). Not used.
655   - * @param[in] phi Azimuth angle in camera system (rad). Not used.
656   - * @param[in] zenith Zenith angle in Earth system (rad). Not used.
657   - * @param[in] azimuth Azimuth angle in Earth system (rad). Not used.
658   - * @param[in] etrue Use true energy (true/false). Not used.
659   - *
660   - * Calculate the radius from the center that contains 'fraction' percent
661   - * of the events. fraction * 100. = Containment % .
662   - ***************************************************************************/
663   -double GCTAPsfKing::containment_radius(const double& fraction,
664   - const double& logE,
665   - const double& theta,
666   - const double& phi,
667   - const double& zenith,
668   - const double& azimuth,
669   - const bool& etrue) const
670   -{
671   - // Update the parameter cache
672   - update(logE, theta);
673   -
674   - // Initialize Radius value
675   - double radius = 0.0 ;
676   -
677   - // use analytic calculation
678   - double arg = std::pow( 1 - fraction, 1/(1-m_par_gamma) ) ;
679   - radius = m_par_sigma * std::pow( 2.0 * m_par_gamma * ( arg - 1.0 ) , 0.5 ) ;
680   -
681   - // Return radius containing fraction of events
682   - return radius;
683   -}
684   -
685   -
... ...
inst/cta/src/GCTAPsfPerfTable.cpp
  View file @2ac11cb
1 1 /***************************************************************************
2 2 * GCTAPsfPerfTable.cpp - CTA performance table PSF class *
3 3 * ----------------------------------------------------------------------- *
4   - * copyright (C) 2012-2013 by Juergen Knoedlseder *
  4 + * copyright (C) 2012-2015 by Juergen Knoedlseder *
5 5 * ----------------------------------------------------------------------- *
6 6 * *
7 7 * This program is free software: you can redistribute it and/or modify *
... ... @@ -37,6 +37,8 @@
37 37  
38 38 /* __ Method name definitions ____________________________________________ */
39 39 #define G_LOAD "GCTAPsfPerfTable::load(std::string&)"
  40 +#define G_CONTAINMENT_RADIUS "GCTAPsfPerfTable::containment_radius(double&,"\
  41 + " double&, double&, double&, double&, double&, bool&)"
40 42  
41 43 /* __ Macros _____________________________________________________________ */
42 44  
... ... @@ -382,6 +384,72 @@ double GCTAPsfPerfTable::delta_max(const double&amp; logE,
382 384  
383 385  
384 386 /***********************************************************************//**
  387 + * @brief Return the radius that contains a fraction of the events (radians)
  388 + *
  389 + * @param[in] fraction of events (0.0-1.0)
  390 + * @param[in] logE Log10 of the true photon energy (TeV).
  391 + * @param[in] theta Offset angle in camera system (rad). Not used.
  392 + * @param[in] phi Azimuth angle in camera system (rad). Not used.
  393 + * @param[in] zenith Zenith angle in Earth system (rad). Not used.
  394 + * @param[in] azimuth Azimuth angle in Earth system (rad). Not used.
  395 + * @param[in] etrue Use true energy (true/false). Not used.
  396 + * @return Containment radius (radians).
  397 + *
  398 + * @exception GException::invalid_argument
  399 + * Invalid fraction specified.
  400 + *
  401 + * Evaluates:
  402 + *
  403 + * \f[
  404 + *
  405 + * radius = \sqrt{ \frac{\ln{\left( \frac{ fraction * m\_par\_width}
  406 + * {\pi*m\_par\_scale} +1\right)}}{m\_par\_width} }
  407 + *
  408 + * \f]
  409 + *
  410 + * which is derived by integrating
  411 + *
  412 + * \f[
  413 + *
  414 + * fraction = \int_{0}^{2\pi}\int_{0}^{radius}r *
  415 + * e^{ m\_par\_width * r^{2}}dr d\phi
  416 + *
  417 + * \f]
  418 + *
  419 + * and solving for radius.
  420 + *
  421 + * Calculate the radius from the center that contains 'fraction' percent
  422 + * of the events. fraction * 100. = Containment % .
  423 + ***************************************************************************/
  424 +double GCTAPsfPerfTable::containment_radius(const double& fraction,
  425 + const double& logE,
  426 + const double& theta,
  427 + const double& phi,
  428 + const double& zenith,
  429 + const double& azimuth,
  430 + const bool& etrue) const
  431 +{
  432 + // Check input argument
  433 + if (fraction <= 0.0 || fraction >= 1.0) {
  434 + std::string message = "Containment fraction "+
  435 + gammalib::str(fraction)+" must be between " +
  436 + "0.0 and 1.0, not inclusive.";
  437 + throw GException::invalid_argument(G_CONTAINMENT_RADIUS, message);
  438 + }
  439 +
  440 + // Update the parameter cache
  441 + update(logE);
  442 +
  443 + // Compute radius
  444 + double arg = fraction * m_par_width / ( gammalib::pi * m_par_scale);
  445 + double radius = std::sqrt(std::log(arg + 1.0) / m_par_width);
  446 +
  447 + // Return containement radius
  448 + return radius;
  449 +}
  450 +
  451 +
  452 +/***********************************************************************//**
385 453 * @brief Print point spread function information
386 454 *
387 455 * @param[in] chatter Chattiness (defaults to NORMAL).
... ... @@ -507,57 +575,3 @@ void GCTAPsfPerfTable::update(const double&amp; logE) const
507 575 // Return
508 576 return;
509 577 }
510   -
511   -
512   -/***********************************************************************//**
513   - * @brief Return the radius that contains a fraction of the events (radians)
514   - *
515   - * @param[in] fraction of events (0.0-1.0)
516   - * @param[in] logE Log10 of the true photon energy (TeV).
517   - * @param[in] theta Offset angle in camera system (rad). Not used.
518   - * @param[in] phi Azimuth angle in camera system (rad). Not used.
519   - * @param[in] zenith Zenith angle in Earth system (rad). Not used.
520   - * @param[in] azimuth Azimuth angle in Earth system (rad). Not used.
521   - * @param[in] etrue Use true energy (true/false). Not used.
522   - *
523   - * Evaluates:
524   - *
525   - * \f[
526   - *
527   - * radius = \sqrt{ \frac{\ln{\left( \frac{ fraction * m\_par\_width}
528   - * {\pi*m\_par\_scale} +1\right)}}{m\_par\_width} }
529   - *
530   - * \f]
531   - *
532   - * which is derived by integrating
533   - *
534   - * \f[
535   - *
536   - * fraction = \int_{0}^{2\pi}\int_{0}^{radius}r *
537   - * e^{ m\_par\_width * r^{2}}dr d\phi
538   - *
539   - * \f]
540   - *
541   - * and solving for radius.
542   - *
543   - * Calculate the radius from the center that contains 'fraction' percent
544   - * of the events. fraction * 100. = Containment % .
545   - ***************************************************************************/
546   -double GCTAPsfPerfTable::containment_radius(const double& fraction,
547   - const double& logE,
548   - const double& theta,
549   - const double& phi,
550   - const double& zenith,
551   - const double& azimuth,
552   - const bool& etrue) const
553   -{
554   - // Update the parameter cache
555   - update(logE);
556   -
557   - // Compute radius
558   - double arg = fraction * m_par_width / ( gammalib::pi * m_par_scale ) ;
559   - double radius = std::sqrt( std::log( arg + 1 ) / m_par_width ) ;
560   -
561   - // Return maximum PSF radius
562   - return radius;
563   -}
... ...
inst/cta/src/GCTAPsfVector.cpp
  View file @2ac11cb
1 1 /***************************************************************************
2 2 * GCTAPsfVector.cpp - CTA point spread function vector class *
3 3 * ----------------------------------------------------------------------- *
4   - * copyright (C) 2012-2014 by Juergen Knoedlseder *
  4 + * copyright (C) 2012-2015 by Juergen Knoedlseder *
5 5 * ----------------------------------------------------------------------- *
6 6 * *
7 7 * This program is free software: you can redistribute it and/or modify *
... ... @@ -37,6 +37,8 @@
37 37 #include "GCTAException.hpp"
38 38  
39 39 /* __ Method name definitions ____________________________________________ */
  40 +#define G_CONTAINMENT_RADIUS "GCTAPsfVector::containment_radius(double&,"\
  41 + " double&, double&, double&, double&, double&, bool&)"
40 42  
41 43 /* __ Macros _____________________________________________________________ */
42 44  
... ... @@ -428,6 +430,72 @@ double GCTAPsfVector::delta_max(const double&amp; logE,
428 430  
429 431  
430 432 /***********************************************************************//**
  433 + * @brief Return the radius that contains a fraction of the events (radians)
  434 + *
  435 + * @param[in] fraction of events (0.0-1.0)
  436 + * @param[in] logE Log10 of the true photon energy (TeV).
  437 + * @param[in] theta Offset angle in camera system (rad). Not used.
  438 + * @param[in] phi Azimuth angle in camera system (rad). Not used.
  439 + * @param[in] zenith Zenith angle in Earth system (rad). Not used.
  440 + * @param[in] azimuth Azimuth angle in Earth system (rad). Not used.
  441 + * @param[in] etrue Use true energy (true/false). Not used.
  442 + * @return Containment radius (radians).
  443 + *
  444 + * @exception GException::invalid_argument
  445 + * Invalid fraction specified.
  446 + *
  447 + * Evaluates:
  448 + *
  449 + * \f[
  450 + *
  451 + * radius = \sqrt{ \frac{\ln{\left( \frac{ fraction * m\_par\_width}
  452 + * {\pi*m\_par\_scale} +1\right)}}{m\_par\_width} }
  453 + *
  454 + * \f]
  455 + *
  456 + * which is derived by integrating
  457 + *
  458 + * \f[
  459 + *
  460 + * fraction = \int_{0}^{2\pi}\int_{0}^{radius}r *
  461 + * e^{ m\_par\_width * r^{2}}dr d\phi
  462 + *
  463 + * \f]
  464 + *
  465 + * and solving for radius.
  466 + *
  467 + * Calculate the radius from the center that contains 'fraction' percent
  468 + * of the events. fraction * 100. = Containment % .
  469 + ***************************************************************************/
  470 +double GCTAPsfVector::containment_radius(const double& fraction,
  471 + const double& logE,
  472 + const double& theta,
  473 + const double& phi,
  474 + const double& zenith,
  475 + const double& azimuth,
  476 + const bool& etrue) const
  477 +{
  478 + // Check input argument
  479 + if (fraction <= 0.0 || fraction >= 1.0) {
  480 + std::string message = "Containment fraction "+
  481 + gammalib::str(fraction)+" must be between " +
  482 + "0.0 and 1.0, not inclusive.";
  483 + throw GException::invalid_argument(G_CONTAINMENT_RADIUS, message);
  484 + }
  485 +
  486 + // Update the parameter cache
  487 + update(logE);
  488 +
  489 + // Compute radius for gaussian
  490 + double arg = fraction * m_par_width / (gammalib::pi * m_par_scale);
  491 + double radius = std::sqrt(std::log(arg + 1.0) / m_par_width);
  492 +
  493 + // Return containement radius
  494 + return radius;
  495 +}
  496 +
  497 +
  498 +/***********************************************************************//**
431 499 * @brief Print point spread function information
432 500 *
433 501 * @param[in] chatter Chattiness (defaults to NORMAL).
... ... @@ -551,58 +619,3 @@ void GCTAPsfVector::update(const double&amp; logE) const
551 619 // Return
552 620 return;
553 621 }
554   -
555   -
556   -/***********************************************************************//**
557   - * @brief Return the radius that contains a fraction of the events (radians)
558   - *
559   - * @param[in] fraction of events (0.0-1.0)
560   - * @param[in] logE Log10 of the true photon energy (TeV).
561   - * @param[in] theta Offset angle in camera system (rad). Not used.
562   - * @param[in] phi Azimuth angle in camera system (rad). Not used.
563   - * @param[in] zenith Zenith angle in Earth system (rad). Not used.
564   - * @param[in] azimuth Azimuth angle in Earth system (rad). Not used.
565   - * @param[in] etrue Use true energy (true/false). Not used.
566   - *
567   - * Evaluates:
568   - *
569   - * \f[
570   - *
571   - * radius = \sqrt{ \frac{\ln{\left( \frac{ fraction * m\_par\_width}
572   - * {\pi*m\_par\_scale} +1\right)}}{m\_par\_width} }
573   - *
574   - * \f]
575   - *
576   - * which is derived by integrating
577   - *
578   - * \f[
579   - *
580   - * fraction = \int_{0}^{2\pi}\int_{0}^{radius}r *
581   - * e^{ m\_par\_width * r^{2}}dr d\phi
582   - *
583   - * \f]
584   - *
585   - * and solving for radius.
586   - *
587   - * Calculate the radius from the center that contains 'fraction' percent
588   - * of the events. fraction * 100. = Containment % .
589   - ***************************************************************************/
590   -double GCTAPsfVector::containment_radius(const double& fraction,
591   - const double& logE,
592   - const double& theta,
593   - const double& phi,
594   - const double& zenith,
595   - const double& azimuth,
596   - const bool& etrue) const
597   -{
598   - // Update the parameter cache
599   - update(logE);
600   -
601   - // Compute radius
602   - // for gaussian
603   - double arg = fraction * m_par_width / ( gammalib::pi * m_par_scale ) ;
604   - double radius = std::sqrt( std::log( arg + 1 ) / m_par_width ) ;
605   -
606   - // Return maximum PSF radius
607   - return radius;
608   -}
... ...
inst/cta/test/test_CTA.cpp
  View file @2ac11cb
... ... @@ -391,10 +391,10 @@ void TestGCTAResponse::test_response_psf(void)
391 391 }
392 392 test_value(sum, 1.0, 0.001, "PSF integration for "+eng.print());
393 393  
394   - // test containment radii
395   - for ( double f = 0.01 ; f < 1.0 ; f += 0.07 ) {
396   - double radius = rsp.psf()->containment_radius( f, eng.log10TeV() ) ;
397   - test_assert( radius>0.0 && radius < 1e6 , "PSF containment radius for fraction="+gammalib::str(f) ) ;
  394 + // Test containment radii
  395 + for (double f = 0.01; f < 1.0; f += 0.07) {
  396 + double radius = rsp.psf()->containment_radius(f, eng.log10TeV());
  397 + test_assert(radius > 0.0 && radius < 1e6, "PSF containment radius for fraction="+gammalib::str(f));
398 398 }
399 399 }
400 400  
... ...