How to calculate event weights for offset ranges using diffuse gammas
Closed this issue · 17 comments
This error appears when using the script calculate_eventdisplay_irfs.py with input files that contain events simulated from different offsets, at line 95:
p["events"]["weight"] = calculate_event_weights(
p["events"]["true_energy"], p["target_spectrum"], p["simulated_spectrum"]
)
The error is: UnitConversionError: 'm2 sr / cm2' (solid angle) and '' (dimensionless) are not convertible
calculate_event_weights tries to divide the target spectrum by the simulated spectrum expecting to get a dimensionless result, but these spectrums have different units in the case of the gammas.
The target_spectrum for the gammas is set to CRAB_HEGRA (from pyirf/spectral.py line 312):
CRAB_HEGRA = PowerLaw(
normalization=2.83e-11 / (u.TeV * u.cm ** 2 * u.s), index=-2.62, e_ref=1 * u.TeV,
)
While the simulated spectrum is gotten at calculate_eventdisplay_irfs.py line 94:
p["simulated_spectrum"] = PowerLaw.from_simulation(p["simulation_info"], T_OBS)
The function 'from_simulation' is defined in pyirf/spectral.py and it sets different units depending on whether viewcone is greater than 0 or not. In this case, as the MCs are from different offsets, the value of viewcone is greater than 0, so the units get an extra u.sr, which produces the error.
Hi,
You are trying to calculate weights for diffuse events given a point-source spectrum. This doesn't make sense and the unit system thus complained correctly.
The Crab Nebula Powerlaw is a point source spectrum (no per steradian).
If you want to have a diffuse flux weighted to something physical, you need a model for the diffuse flux
@GernotMaier Could you please explain how do you perform the weighting of the events for diffuse gammas? How do you include the sr unit into the flux?
This needs to be implemented within pyirf for calculating full-enclosure IRFs.
I think if the task is to assign Crab-like weights to diffuse gamma-ray simulations between some field of view offset, you need to integrate the diffuse flux you get from the simulation for this and then use the resulting bit for the event weights, like this:
import astropy.units as u
from pyirf.simulations import SimulatedEventsInfo
from pyirf.spectral import CRAB_HEGRA, PowerLaw, POINT_SOURCE_FLUX_UNIT, calculate_event_weights
from pyirf.utils import cone_solid_angle
# calculate sensitivity between 1 and 2 degrees offset from fov center
obstime = 50 * u.hour
theta_low = 1 * u.deg
theta_high = 2 * u.deg
# just an example, in reality, read from simulations
simulated_events = SimulatedEventsInfo(
n_showers=int(1e6),
energy_min=10 * u.GeV,
energy_max=100 * u.TeV,
max_impact=1 * u.km,
spectral_index=-2,
viewcone=10 * u.deg,
)
simulated_spectrum = PowerLaw.from_simulation(simulated_events, obstime=obstime)
# "integrate" the spectrum between theta_low, theta_high
solid_angle = cone_solid_angle(theta_high) - cone_solid_angle(theta_low)
spectrum_ring = PowerLaw(
(simulated_spectrum.normalization * solid_angle).to(POINT_SOURCE_FLUX_UNIT),
simulated_spectrum.index,
simulated_spectrum.e_ref,
)
print(calculate_event_weights(5 * u.TeV, CRAB_HEGRA, spectrum_ring))Hi @JBernete,
the underling issue here is that, even if technically the IRF implementation in pyirf is already able to write IRFs which do not depend on the offset angle, we never actually tested pyirf on diffuse simulation as we decided to validate the project against EventDisplay from the easiest case first (best point-source flux sensitivity).
I hope that this issue will now trigger this missing validation and you are of course welcome to contribute to the project accordingly.
@HealthyPear the issue here has nothing to do with IRFs, it is about how to use diffuse simulations to calculate a point source sensitivity.
the issue here has nothing to do with IRFs
@maxnoe I feel I disagree with this. Yes, technically it has nothing to do with IRFs. But effectively, the first thing IRF producers look at after a production, is sensitivity. And definitely not only point-source sensitivity: the main objective for IRF production is for science tools to use them, and science tools require full-enclosure IRFs for most cases.
So I would say ensuring we are able to validate and cross-check eventDisplay full-enclosure IRF production, of course using diffuse gammas, should be the next milestone of pyirf, and that (to my knowledge) has not been done yet.
@TarekHC fully agree on the general "diffuse sensitivity / irf part". I meant this specific issue about calculating weights for diffuse events.
Yes, technically it has nothing to do with IRFs
Exactly, please, let's not mix technical issues that might be solved quickly with broad, general discussion.
To be fair, I never said the problem are the IRFs...I said that the problem is that we never validated pyirf against EventDisplay diffuse case before.
Let's discuss the general picture in a dedicated issue and solve the technical problem at hand here.
We mix to many things in a single issue, it will make it hard to follow and later find the discussion.
I opened #195 for the general discussion. Let's focus here on the event weights and wait for @JBernete to comment if the example I showed produced the expected results or @GernotMaier to comment how this is handled in EventDisplay.
On the technical side, it might be nice to implement what I wrote above as a PowerLaw.integrate_solid_angle(theta_low, theta_high) -> PowerLaw method.
I am not sure I understand the discussion in all details. IRFs are calculated assuming point-like sources at different offsets and the calculation is identical to the on-source case. The only difference is that the approximation of a point-like source at a given offset by using MC events from a offset interval (approximation: angular differences are calculated relative to the simulated direction, not relative to the source position), otherwise the calculation is exactly the same (counting the number of events after cuts vs the number of simulated events).
The solid angle is relevant only for the proton/electron weighting.
The solid angle is relevant only for the proton/electron weighting
So you use histograms of the simulated CORSIKA events in energy/offset for the weighting?
In that case, we still need the solid angle for the gammas, since we (at least for now due to the problem with the exported histograms) use the production spectrum of CORSIKA and the total number of events to calculate that.
Exactly, I fill histograms as function of offset and energy with the simulated events (with the disadvantage that this needs to be done at a very early stage of the analysis and the binning is kind of fixed then).
And yes, using the production spectrum + viewcone definition from CORSIKA requires taking the solid angle normalisation into account.
You could probably get around the first issue by interpolation the histogram assuming a powerlaw inside the bins... which we might need to do anyway here, since we will have only the prebinned histograms of simtelarray in the ctapipe output
I opened #195 for the general discussion. Let's focus here on the event weights and wait for @JBernete to comment if the example I showed produced the expected results or @GernotMaier to comment how this is handled in EventDisplay.
On the technical side, it might be nice to implement what I wrote above as a
PowerLaw.integrate_solid_angle(theta_low, theta_high) -> PowerLawmethod.
Yes, your example is producing the expected results. Thank you. I also think it would be a good idea to implement it as a method for PowerLaw.