|Assigned To:||Mayer Michael||% Done:|
To be able to reproduce results of current IACTs, we need to implement an elliptical Gaussian model.
#1 Updated by Deil Christoph over 7 years ago
Is it necessary to duplicate the radial and elliptical models in Gammalib?
E.g. Sherpa only has a Gaussian model, it has elongation parameters, the default is that they are frozen, i.e. the user gets a radial Gaussian model:
This page also shows the formulas to implement the elongation in a generic way that could be re-used by all elongated models.
But maybe a single efficient implementation is no longer possible because some analytical integrals are needed for unbinned analysis?
#2 Updated by Knödlseder Jürgen over 7 years ago
There could be eventually some re-organising in the future of these models, but for the time being we have separated radial from elliptical models. Doing now a Gaussian in addition to a disk model is a no brainer, it should even work much better than the disk as it has no sharp edges.
I also see one difference between radial and elliptical models: for radial models I’m sure they are done correctly on the sphere (hence they could be used for example for a Fermi-LAT analysis to describe Loop I). For elliptical models I’m not really sure in what system the computation is done, I think it’s cartesian.
#3 Updated by Mayer Michael over 7 years ago
- File pulldist.png added
- Status changed from New to Feedback
- Assigned To set to Mayer Michael
- Target version set to 1.0.0
- % Done changed from 0 to 90
I have implemented the new class
GModelSpatiallEllipticalGauss, which was actually just copying the majority of the code from the elliptical disk model. I still have to update some parts of the documentation I guess (which is quite hard since the formula is quite lengthy).
Nevertheless, the code seems to work properly. However, I encountered rather long computation times when I run
cspull. By constraining the parameter boundaries closer to the real values some speed could be gained but an fit of a “standard” CTA observation still takes about 20 secs on my MacBook.
Therefore, I reduced the observation time in
cspull to 900sec (instead of 1800). I quickly ran 200 trials and created some pull distribution for the parameters (which are attached).
I encountered one issue, we might want to think about:
In the base class
GModelSpatialElliptical, the parameter names of the axes are set to “MinorRadius” and “MajorRadius”. For the Gaussian, I’d rather go for “MinorAxis” and “MajorAxis” since “radius” is not the proper name when it comes to a Gaussian. Should we therefore reorganise the elliptical classes to move the two size parameters to the derived classes, or should
GModelSpatialEllipticalGauss::init_members() just rename the parameters?
I furthermore created a test case, which passes all tests (analogous to the disk model).
Everything is available on branch 1412-implement-GModelSpatialEllipticalGauss
#4 Updated by Deil Christoph over 7 years ago
Sherpa parametrises the elongated models like 2D Gaussian with width and ellipticity (see http://cxc.harvard.edu/sherpa/ahelp/gauss2d.html).
Astropy uses semi-major and semi-minor axis like you did here.
Does one of you know which one is more stable for morphology fitting?
I remember when I tried to make a HESS source catalog of elongated Gaussians with Sherpa getting the fits to converge in complex regions was a big problem even with good starting values for the optimizer and I eventually gave up.
Do you have a test that shows that the position angle is fit robustly near the theta = 0 deg = 360 deg ?
In any case, I’m +1 to put in what you have now in devel and maybe encourage people to try out elongated models in the mailing list, so that potential issues are noticed before the 1.0 release.
#5 Updated by Knödlseder Jürgen over 7 years ago
We may indeed change the parameters, one could see what is more stable.
There is in fact a fundamental problem with major and minor axis, which is that you cannot guarantee that at the end the major axis is in fact the longer of both axes (a fit may drive it to a smaller value than the minor axis).
With ellipticity one can constrain the parameter in a way that the major axis is always the longer axis.
#6 Updated by Deil Christoph over 7 years ago
For reference, here’s the parametrisation Astropy uses:
The also have the option to initialise the parameters from a matrix, but that’s not the parametrisation used in the fit.
If I had to guess and decide on a parametrisation without a study (which is not trivial to do to cover all typical CTA use cases for fitting elongated models to see if it’s stable), I’d pick the Sherpa (width, elongation) parametrisation.
#7 Updated by Mayer Michael over 7 years ago
Thanks for the links Christoph. I used exactly the same formula as astropy.
Changing the parameters to width and ellipticity should probably already happen in the
GModelSpatialElliptical base class, right? Then we would have to change the elliptical disk model accordingly.
If we agree that we want to go that way, I can try to apply these changes. What do you think?
Jürgen, should that happen in the branch for this issue, or is this something which deserves another issue/change request?
#8 Updated by Deil Christoph over 7 years ago
Since probably none of us will have time for a detailed study, I’ve asked for advice here:
I’m +1 to merge whatever you have now in devel since it’s not clear if we should change to another parametrisation ... increases the chances of me and other trying it out on a simulated CTA GPS or the real HESS GPS data.
#9 Updated by Deil Christoph over 7 years ago
I’ve discussed this with Axel a bit.
Probably both the Astropy and Sherpa parametrisation is unstable for roughly circular sources,
because the position angle jumps when one of the two width parameters gets larger than the other.
This might have been the problem when I tried to make a HESS survey catalog of elongated Gaussians,
and in any case for CTA Galactic analysis we want a tool gives robustly converging elongated model fits for all sources, even the ones that are roughly circular.
One parametrisation that might be stable is CXX, CXY and CYY from SExtractor, see the section 10.6.1 in the manual:
It might be worthwile to use this internally or to implement this as an alternative parametrisation.
This is similar to spectral models where in some parametrisations problems occur, e.g. in the HESS software we have 2 power-law parametrisations and 3 exponential-cutoff power-law parametrisations implemented.
#12 Updated by Deil Christoph over 7 years ago
+1 to call the axes “A”, “B” or “1”, “2” instead of “major” and “minor” and then explain via documentation that the user has to look afterwards themselves which one is the major and minor one.
PS: the latest version of the branch can be seen here:
#13 Updated by Knödlseder Jürgen over 7 years ago
Before making any decision I would like to analyse the real stability of the actual code (for example to see how an ellipse is fitted where the PA starting value differs by 45 degrees, etc.)
Alternatively, I’d like to try replacing
minor_axis^2 = major_axis^2 (1 - eccentricity^2). This allows to constrain the
eccentricity to the interval [0,1], making sure that the
major_axis is always the major axis. I’d like to compare the stability of that approach to the actual implementation.