GammaLib

Comparing the psf in different classes, obserations, telescopes, and code formats can difficult, but the psf’s containment radius is a pretty universal way to do this. However, it seems theres no native way in gammalib/ctools to calculate this, making comparisons and error checking difficult.

So, I propose an additional virtual function for GCTAPsf (and its inheriting classes), where a user gives a containment fraction(0.0-1.0), and the point in the energy/theta/phi/zenith/azimuth parameter space, and the containment radius containing that fraction of events is returned, in radians. I imagine the function prototype in GCTAPsf.hpp would look something like:

    virtual double      containment_radius( const double& fraction,
                                  const double& logE,
                                  const double& theta = 0.0,
                                  const double& phi = 0.0,
                                  const double& zenith = 0.0,
                                  const double& azimuth = 0.0,
                                  const bool&   etrue = true) const = 0;

I have written a python function that calculates the containment radius from the existing GCTAPsf::operator(), but it has to take a bunch of samples(slow), and then linearly interpolate between the two closest ones(inaccurate). Theres probably a better algorithm, or it may be directly computable from the psf function parameters, but the code below may still be useful to whomever is adding these functions.

import gammalib, numpy

def containment_radius( run, en, th, contain=0.5, npts=100 ) :
  """Calculate the containment radius that contains 'contain' fraction of the events.

  Args:
    contain: float, 0.0-1.0, containment fraction to find radius for
    run    : GCTAObservation object, should have psf irf loaded
    en     : float, energy, log10TeV
    th     : float, camera offset, radians
    npts : number of times to sample the psf between 0 and delta_max, higher npts = higher accuracy in the returned value

  Returns:
    float, radians, radius containing 'contain' fraction of events
  """ 
  samps = numpy.zeros(npts)
  dmax  = run.response().psf().delta_max( en, th ) # radians
  delta = dmax / npts # radians

  # loop over each point
  for i in range(npts) :
    d = i * delta # radians
    counts_per_sr = run.response().psf()( d, en, th )   # counts per sr
    if counts_per_sr == 0.0 :
      print 'warning: containment_radius(en=%f, th=%f): 0 counts at radius %.3f deg, probably no psf info' % ( en, th, d * gammalib.rad2deg )
      return 0.0

    # integrate these counts around the circle of radius d
    # counts * perimeter of circle with that radius, 
    # may get less accurate when d is large (>>10deg?)
    samps[i]      = counts_per_sr * 2 * gammalib.pi * d

  # total counts in our psf from 0.0 to dmax
  total_counts = numpy.sum(samps)

  # loop over each sample, check if we cross the containment fraction
  for j in range(npts-1) :
    partial_counts      = numpy.sum( samps[:j+1] )
    partial_counts_next = numpy.sum( samps[:j+2] )
    cfrac      = partial_counts      / total_counts
    cfrac_next = partial_counts_next / total_counts
    radius      =  j    * delta
    radius_next = (j+1) * delta

    # check if we have crossed the containment fraction
    if cfrac < contain and cfrac_next > contain :
      # if so, linearly interploate to find what radius has the target containment fraction
      # starting from http://en.wikipedia.org/wiki/Linear_equation , the Two-point form
      contain_radius = radius + ( (contain - cfrac) * (radius_next - radius) / (cfrac_next - cfrac ) )
      return contain_radius

  raise LookupError('Could not find containment fraction %.2f within radii 0.0 and %.3f (radians)' % (contain, dmax) )