Implement the Earth Horizon cut

[GallegoSav]

This pull request add the Earth Horizon cut feature. This is done by 2 steps :

1. Match the _ez_cart and _min_angle_cos computed by get_earth_occ for each event. Then compute the angle between PsiChi and _ez_cart .
2. Compute the fraction of the Compton cone that is above the Earth Horizon (0 completely below, 1 completely above)

The user can then set a cut value based on this distribution with cut_EarthHorizon().
By using only numpy, this function is pretty fast ( ~3,4 min for 118 M events of the Cosmic photons bg component).

Here an example on how this distribution differs between the Albedo Photons (coming from below the satellite) and the Cosmic Photons (coming from above).

The power of discrimination is not that great but that's not a surprise. This can be of course improve by asking nb of interaction > 2 and applying some cut in Phi, if the signal is a point source , ect...

[read-the-docs-community]

Documentation build overview

📚 cositools-cosipy | 🛠️ Build #33164429 | 📁 Comparing 99009d0 against latest (4dd844d)

  🔍 Preview build  

14 files changed · ± 9 modified · - 5 deleted

± Modified

• genindex.html
• api/index.html
• api/polarization.html
• api/threeml.html
• tutorials/DataIO/DataIO_example.html
• tutorials/polarization/ASAD_method.html
• tutorials/polarization/Stokes_method.html
• tutorials/polarization/maximum_likelihood_method.html
• tutorials/source_injector/Point_source_injector.html

- Deleted

• py-modindex.html
• _modules/index.html
• _modules/cosipy/polarization/conventions.html
• _modules/cosipy/polarization/polarization_angle.html
• _modules/cosipy/polarization/polarization_axis.html

[codecov]

Codecov Report

❌ Patch coverage is 70.34314% with 121 lines in your changes missing coverage. Please review.
✅ Project coverage is 77.44%. Comparing base (e3d4d83) to head (99009d0).

Files with missing lines	Patch %	Lines
cosipy/data_io/EmCDSUnbinnedData.py	58.04%	86 Missing ⚠️
cosipy/event_selection/earth_horizon_selection.py	75.96%	25 Missing ⚠️
cosipy/data_io/UnBinnedData.py	90.74%	5 Missing ⚠️
cosipy/interfaces/data_interface.py	84.61%	4 Missing ⚠️
cosipy/interfaces/event.py	94.73%	1 Missing ⚠️

Files with missing lines	Coverage Δ
cosipy/interfaces/event.py	91.37% <94.73%> (+0.68%)	⬆️
cosipy/interfaces/data_interface.py	66.17% <84.61%> (+3.56%)	⬆️
cosipy/data_io/UnBinnedData.py	92.26% <90.74%> (-0.27%)	⬇️
cosipy/event_selection/earth_horizon_selection.py	75.96% <75.96%> (ø)
cosipy/data_io/EmCDSUnbinnedData.py	52.00% <58.04%> (+6.91%)	⬆️

🚀 New features to boost your workflow:

• ❄️ Test Analytics: Detect flaky tests, report on failures, and find test suite problems.

[israelmcmc]

Thanks @GallegoSav. It's fine for now, but it would be good to make it an EventSelector instead. We plan to refactor the DataIO anyway, see #344

[jdbuhler]

Quick review, mainly for performance issues since you indicated that you were interested in speed on a large set of events:

• it might be preferable to take the SpacecraftHistory as an argument rather than re-opening the file each time, which has a nontrivial cost.
• Angle_PsiChi_Ez() returns two angles, but at line 946/947, you call the first return value costheta?
• Line 990: what is the "TimeTags" element of the dict before you sort it? You're converting to an array, which makes a copy
• Lines 1009-1011 recompute the same arrays of cosines and sines multiple times
• A quick google suggests that NumPy's einsum, unlike dot and tensordot, does not call OpenBLAS to accelerate multiplication or use multiple threads. Is there a way to express line 1031 with tensordot?
• arccos is expensive. The angles returned from Angle_PsiChi_Ez() are immediately passed into calculate_sky_fraction(), where they are reduced to their sines/cosines immediately. Can conversion back to angles be avoided?
• by the same token, can the arccos to compute f_sky be avoided if you instead take the cosine of cutvalue?

[jdbuhler]

(BTW, if the lack of acceleration for einsum is a real thing, we should be worrying about it elsewhere in cosipy as well. It comes up in relative coordinates conversion, location and attitude interpolation in SpacecraftHistory, and whatever is happening in optimized_unbinned_folding.py)

[jdbuhler]

OK, after some further investigation, it's not clear that the einsums noted above are bottlenecks, or that they can easily be replaced by dot/tensordot calls. So maybe never mind on that point.

(More specifically, the operation in this code is equivalent to np.vecdot(), but the latter, despite being advertised as using CBLAS under the hood, refuses to use more than one core no matter how large I make either the number of vectors or the vector length. vecdot() is actually somewhat slower on my machine than the corresponding einsum() call, which is appalling. I can get the expected performance with multiple cores if I explicitly rewrite the vecdot() operation with Numba, but I'm not sure it's worth it here.)

[ckarwin]

@GallegoSav I'm happy to start looking at this next week if somebody hasn't approved it yet.

[GallegoSav]

@israelmcmc @jdbuhler Thanks for your comments ! I will try to make it an EventSelector .

[GallegoSav]

@israelmcmc @jdbuhler I implemented a EventSelectorInterface of the Earth Horizon cut.

However, I get small differences in the distributions of fsky between the 2 methods. It could come from :

1. With the interfaces I am now matching the galactic pointing for each events using the spacecraft attitude, where before I was directly using PsiChi galactic .
2. I need to convert the ori.obstime (unix) into jullian date format in order to compare with the time format used in the event interfaces.

[israelmcmc]

Interesting. Assuming there not an error somewhere, I think the largest difference will come from the PsiChi local to Galactic conversion. It would be good to distribution of the difference between these 2 ways to convert PsiChi. I think @ckarwin looked into this at some point. @GallegoSav are you using the 1s cadence ori file or larger?

This and issues that have come up have convince me that I would be move convenient and efficient if PsiChi were available in both coordinates systems for every event, instead of computing it on the fly every time. This is true for UnbinnedData, but not yet for the interfaces. We would need an equivalent to ComptonDataSpaceInSCFrameEventInterface  --e.g. ComptonDataSpaceInGalacticEventInterface and make a derived interface that inherits from both.

[GallegoSav]

Interesting. Assuming there not an error somewhere, I think the largest difference will come from the PsiChi local to Galactic conversion. It would be good to distribution of the difference between these 2 ways to convert PsiChi. I think @ckarwin looked into this at some point. @GallegoSav are you using the 1s cadence ori file or larger?

This and issues that have come up have convince me that I would be move convenient and efficient if PsiChi were available in both coordinates systems for every event, instead of computing it on the fly every time. This is true for UnbinnedData, but not yet for the interfaces. We would need an equivalent to ComptonDataSpaceInSCFrameEventInterface --e.g. ComptonDataSpaceInGalacticEventInterface and make a derived interface that inherits from both.

@israelmcmc I tried the 1s ori file but it did not change the result. Yes, my first idea was to make this ComptonDataSpaceInGalacticEventInterface class actually and use directly the PsiChi galactic from the event like the unbinned class method.

I added new interfaces for this pull request, where you can have local and galactic PsiChi using TimeTagEmCDSEventDataInSCAndGalFrameFromDC3Fits.

It reduced the differences but that is still not perfect. Maybe the conversion from unix time to julian date is making the last differences ?