GammaLib

+1 to having this in gammalib.

Manuel, you are aware of the HESS internal note on the PSF and morphology fitting, right?
https://bitbucket.org/hess_software/hess_psf/downloads/HESS_PSF_Morphology_2013-02-07.pdf
https://bitbucket.org/hess_software/hess_psf/wiki/PSF_Models
(this is HESS internal because it contains confidential information)

I looked into the one-parameter Student T distribution Kora used in her thesis and found out that it is actually a special case of the two-parameter King profile that Fermi uses, and that the King profile is normalized to integrate to one, which is important for the PSF.

I think we should implement the King PSF (or make it usable for HESS if it is already there) and then the Student T will just be a special case.

See the attached IPython notebook that I used to learn some basic things about PSFs.

Here is the formula I implemented for KING a week ago (it looks different because it is the probability density in the theta square histogram that we usually fit):
http://www.mpi-hd.mpg.de/hfm/HESS/intern/Software/SASH-Docs/PSFInfo_8C_source.html

00614   else if( fitfunction == "TRIPLEEXP" ){
00615 
00616     func =  new TF1( fname.c_str(),
00617              "[0]*(exp(-x/(2*[1]*[1]))+[2]*exp(-x/(2*[3]*[3]))+[4]*exp(-x/(2*[5]*[5])))",
00618              0.000, 0.05 );
00619 
00620     double m = psf->GetMaximum();
00621     func->SetParameters( m, 0.01, 0.1, 0.1, 0.1, 0.2 );
00622     func->SetParLimits(0,0.001,10000.);
00623     func->SetParLimits(1,0.0,1.0);
00624     func->SetParLimits(2,0.001,10000.);
00625     func->SetParLimits(3,0.0,1.0);
00626     func->SetParLimits(4,0.001,10000.);
00627     func->SetParLimits(5,0.0,1.0);
00628     func->SetParNames("scale", "#sigma_{1}", "A_{2}", "#sigma_{2}", "A_{3}", "#sigma_{3}");
00629   }
00630   else if( fitfunction == "KING" ){
00631 
00632     func =  new TF1( fname.c_str(),
00633              "[0] / (2 * [1] * [1]) * (1 - 1 / [2]) * pow(1 + 1 / (2 * [2]) * x  / ([1] * [1]), -[2])",
00634              0.000, 0.05 );
00635 
00636     func->SetParNames("scale", "sigma", "gamma");
00637     func->SetParameter("scale", 1);
00638     func->SetParLimits(0, 0.1, 10.);
00639     func->SetParameter("sigma", 0.05);
00640     func->SetParLimits(1, 0.0, 1.0);
00641     func->SetParameter("gamma", 1.5);
00642     func->SetParLimits(2, 1, 10.);
00643   }

Attached find a short description of the King profile from the internal note.

Manuel and Jürgen: do you know if there are analytical formulae for containment fraction (for a given radius) and containment radius (for a given fraction)
for the triple-Gauss and King profile? That would be useful to have for me. (and I guess also for gammalib?)

For a single Gauss there are simple analytical formulae for containment fraction and radius (see attached screenshot gauss_containment.png).
No numerical integration needed to compute those.

For the multi-Gauss I didn’t have time to think about it yet, if one of you figures it out (or knows of a reason why there shouldn’t be a formula), please let me know.

First, Christoph I think you mean Michael instead of Manuel, right?
Of course I knew about the PSF internal note... In fact, it was Kora’s suggestion to implement the Student’s T function when she saw my residuals. I absolutely agree, that since this function is just a special case of the King Profile, we should implement this instead. Hence, I changed the topic of this thread.
The King Profile can be integrated analytically. But I am still struggling to derive the parameter sigma for a given R68.
In my opinion, the following equation has to be solved to find sigma:
In the HESS lookup tables R68 is calculated without parameterising the PSF. It is done by integrating over Monte Carlo Theta2-distributions for point sources (as far as I know). Therefore, we need to calculate the corresponding sigma to a given R68 (and gamma). So far, I didn’t find an analytical way to do this.
I already started to create the files GCTAKing.hpp and .cpp. I would use GCTAPsf instead of GCTAPsfVector as a base class, if that is fine with you.

Yes, I want to find a formula for sigma(gamma,R68,F) with F being 0.68. I did some simulations for various parameters of gamma, sigma and r68, maybe to find an easy relation between those. In fact, for a fixed gamma, the relation between sigma and r68 couldn’t be easier: It is just a straight line from the origin. The slope of this line again depends on gamma (see plots). The relation between gamma and the slope is a bit more complex but also well fit by a logarithmic function. I found the following relation between the slope S and gamma:

Hence, we could easily parameterise sigma(gamma,R68) at least for the case of F = 0.68 and gamma between 1.2 and 2.5 (I think this might also be the range of gamma in the Fermi PSF?).

I am not sure wether we also need this also dependent on the containment fraction F since R68 is the most commonly used value (maybe we could also do for R95?). By the way: Is F=0.68 exactly, or is the real value more like 0.6827?

Attached a document that provides the analytic formulae: gammalib_maths.pdf.

The linearity between sigma and r68 is seen in Eq. (9) of section 2.1. The slope is given by Eq. (8), which differs from your equation. Below the graphical representation of Eq. (8) to compare with your slope plot. For large gammas the plot seems to fit yours, but for small gammas I find smaller values than you do. Could this be a numerical problem in your estimate, as for small gammas the r68 becomes very large?

For the record, here is a small Python script that checks the analytical formulae: king.py

Just to give a short update: The code is working. I still have to test the mc()-method. If this was done successfully, I will push the branch into the repository. As input format, I followed a GCTAResponseTable using a GFitsTable as input.

Merged into devel