GammaLib

Jürgen, FYI, I am considering switching the HESS Galactic plane survey catalog pipeline from Sherpa to gammalib.
What is missing as far as I see is elongated source models and on-off-likelihood.
I have a student starting next week and hope we’ll finally get around to implementing the things we talked about last fall for gammalib, but it’s not clear if we will.

So +10 on implementing this in gammalib from my side.

Looks good. And I have nothing against being ambitious Once this is in place it will be fantastic.

Only for the Gauss, I propose we add a second s at the end for consistency.

• 2 for the centre position
• 1 for position angle
• (at least) 2 for width along and perpendicular to the position angle (an elliptical shell will be more complicated)

I agree, Gauss and disk are the most (only?) commonly used elongated models, so that’s the place to start.

Before you invest a lot of time: are you sure that it’s better to have separate classes instead of only having elongated Gaussians and disks where the symmetric versions are a special case?
This is how Sherpa and HAP / ParisAnalysis in the HESS software all do it:
http://cxc.harvard.edu/sherpa/ahelp/gauss2d.html

But I guess in the case of gammalib having separate classes would be better?
(I don’t understand how unbinned fitting on the sphere works, so I have no idea what math / computations this implies.)

I absolutely agree to add an s to Gaus
Yes, the elongated Gaussian would have these 5 parameters. Maybe for this one we can learn from HESS and/or Sherpa. I didn’t even think about an elliptcial shell, but this would be very interesting. Indeed, these radial models are special cases of the elliptical ones, which would justify inheritations. But I don’t have a strong opinion on that.

I think it’s better to have separate classes, since the computations will be quite different at the end. Having a special elliptical categorie will allow implementing optimised methods for elliptical models, while the radial models will have their own optimizations. Also, the XML parameter interface would be different for symmetric and elliptical models, and this is easier to handle by specific classes.

But some more thinking is needed for the global model organisation, as the radial models are actually the only possible spatial component of extended models. Probably we can add a branching point here (allowing the spatial component to be of type radial or elliptical), but I need to check the code and interfaces to see whether this needs some modifications.

Concerning the unbinned fitting, the models have to be integrated over a region of interest. The region of interest is a circle, and for radial models I had to work out the intersection of the circular model with the circular ROI (hence the intersection between 2 circles needed to be found). Now we have to work out the intersection of an ellipse with a circle, and then make sure that the integral is only taken over the intersection region.

This sounds like a lot of work to make unbinned fitting with elliptical models work.
For the HESS Galactic plane survey I use ROIs of approx. 10 x 5 deg, i.e. rectangles where I fit groups of overlapping sources.
So this would not be possible with gammalib because only circular ROIs are supported?

In the end, do we need unbinned fitting?
Some reference points: Sherpa only has binned fitting and as I understand it it’s also out of fashion for Fermi, it was never available in the HESS software.

On the other hand, unbinned fitting was initially the cool thing I was excited about and looked at gammalib.
And probably if per-event position errors can be taken into account in the likelihood function somehow, a quite large gain in sensitivity could be obtained.
Maybe not everything (like elliptical models or rectangular ROIs) has to be possible in unbinned mode?

Well, it’s some work, but I won’t qualify this as a lot But I think it’s worth the investment (you have to sweat just once, when it’s done it’s pretty cool, as you say).

In fact, the npred method has even not to be implemented, as a default method exists in the base class, see GResponse::npred_extended. The only reason for overloading this method and for calculating the intersection between the ROI and the model is efficiency.

Concerning the circular ROI, there is one ROI per observation (or run), and since the field of view is basically circular, a circular ROI should be fine. I guess your 10x5 deg regions are after combining the data, but for unbinned, nothing is combined.

I added a couple of subtasks that specify more precisely what needs to be done.

I implemented now the handling of elliptical source models in the CTA response class.

What is missing now is the coding of an elliptical source model to test the code. Michael, are you willing to implement an elliptical disk model? See action #761.

Once we have this model we can test the code and make sure that it works correctly.

Yes, I’ll start.