Add function to compute bias and resolution from energy dispersion #144
Conversation
Codecov Report
@@ Coverage Diff @@
## master #144 +/- ##
==========================================
- Coverage 87.35% 84.95% -2.40%
==========================================
Files 38 38
Lines 1329 1363 +34
==========================================
- Hits 1161 1158 -3
- Misses 168 205 +37
Continue to review full report at Codecov.
|
There was a problem hiding this comment.
for the on-axis case we shouldn't expect big differences
I think doing this from the IRF instead of from the events themselves is more general and the energy benchmarks (meaning res and bias) would get the validation from the IRF indirectly
perhaps (apart from the ED comparison) we could think of dump the old function in favor of this one
|
@HealthyPear I don't really see a reason to dump the old function, both are fine and needed I'd say. |
HealthyPear
added
benchmarking
| bias = np.zeros((n_energy_bins, n_fov_bins)) | ||
| resolution = np.zeros((n_energy_bins, n_fov_bins)) | ||
|
|
||
| for energy_bin in range(n_energy_bins): |
There was a problem hiding this comment.
The only way I have found in Gammapy to avoid the loop is the "brute-force" approach of finding the closest corresponding value: compute the cdf and setup an ND interpolator, evaluate the cdf on a finer grid, with the precision you require and then use np.argmin(cdf_fine - value, axis=) to find the closest value along the given axis. The resulting ND implementation uses a lot numpy broadcasting, see https://github.com/gammapy/gammapy/blob/master/gammapy/irf/psf/core.py#L53, for a similar problem. I'm not sure whether it is more efficient in the end, but the solution is at least ND and does not rely on the hard-coded loops over the other dimensions. Not sure I would recommend it though...
There was a problem hiding this comment.
I just remember @registerrier implemented an ND solution, based on np.nditer() for a similar problem to avoid the hard-coded Python loops here: https://github.com/gammapy/gammapy/blob/master/gammapy/stats/counts_statistic.py#L42 Might be useful as well..
There was a problem hiding this comment.
The advantage of the first solution is that it fully supports numpy broadcasting. So something like this:
containment = np.array([0.68, 80, 90])
energy_true = [1, 3, 10, 30, 100] * u.TeV
offset = [0.5, 1] * u.deg
psf.containment_radius(
containment=containment.reshape((-1, 1, 1)),
energy_true=energy_true.reshape((1, -1, 1)),
offset=offset.reshape((1, 1,- 1))
)Returns the expected broad-casted result. Not sure it's needed here but in general it's nice.
No description provided.