Bug in I3HighestEparticleExtractor

[sevmag]

Describe the bug

This issue lists a few bugs of the I3HighestEparticleExtractor

1. Negative visible length

The visible length as computed for tracks can be negative if both intersection points are negative. To my knowledge, this happens when the particle never intersects the hull. (I'm not too sure about that)

graphnet/src/graphnet/data/extractors/icecube/i3highesteparticleextractor.py

Lines 336 to 338 in f3faff1

	visible_length = intersections.second - max(
	intersections.first, 0
	)

2. Track filtering only for top-level primaries

Currently, filtering for tracks that are in the subtrees of the primaries is handled here:

graphnet/src/graphnet/data/extractors/icecube/i3highesteparticleextractor.py

Lines 186 to 192 in f3faff1

	MMCTrackList = [
	track
	for track in MMCTrackList
	if frame[self.mctree].get_primary(track.GetI3Particle())
	in primaries
	]
	MMCTrackList = simclasses.I3MMCTrackList(MMCTrackList)

The primaries are coming from the get_primary function, which can also return particles that are not in the head of the MCTree.

the get_primaries funciton:

graphnet/src/graphnet/data/extractors/icecube/i3extractor.py

Line 166 in f3faff1

def get_primaries(

In the scenario that the returned primaries are not in the head of the MCTree, the filtering logic will always filter out the tracks even if they are descendants of particles in the primaries list. This is because frame[self.mctree].get_primary() will only give particles at the head level.

same issue occurs here:

graphnet/src/graphnet/data/extractors/icecube/i3highesteparticleextractor.py

Lines 399 to 414 in f3faff1

	e_p = np.array(
	[
	np.array([p.energy, p])
	for p in particles
	if (
	(
	(p.energy > min_e)
	& (
	frame[self.mctree].get_primary(p.id)
	in primaries
	)
	)
	& (not p.is_track)
	)
	]
	).T

This problem leads to the energy proxies calculated by this extractor to be 0 or inaccurate if the particles in the primaries are not in the top level of the MCTree.

3. Wrong return in get_bundle_HEP

This returns a single particle

graphnet/src/graphnet/data/extractors/icecube/i3highesteparticleextractor.py

Line 232 in f3faff1

return dataclasses.I3Particle()

However, the function is expected to return a Tuple as seen here:

graphnet/src/graphnet/data/extractors/icecube/i3highesteparticleextractor.py

Line 672 in f3faff1

bundle, length_mask, no_intersect = self.get_bundle_HEP(

I have fixed these bugs already. Positing this here for completeness. @Aske-Rosted, you might want to take a look again.

[sevmag]

Update:

There is another issue with the HEP extractor for bundles:

graphnet/src/graphnet/data/extractors/icecube/i3highesteparticleextractor.py

Lines 695 to 725 in f3faff1

	for track, MGtrack in zip(MMCTrackList, MuonGun_tracks):
	intersections = self.hull.surface.intersection(
	MGtrack.pos, MGtrack.dir
	)
	track_energy = MGtrack.get_energy(intersections.first)
	EonEntrance += track_energy
	closest_pos.append(
	np.array(
	[
	track.GetXc(),
	track.GetYc(),
	track.GetZc(),
	]
	)
	* track_energy
	)
	if closest_time is None:
	closest_time = track.GetTc()
	elif closest_time < track.GetTc():
	closest_time = track.GetTc()
	if intersections.second > 0:
	visible_length = max(
	visible_length, intersections.second - intersections.first
	)
	if MGtrack.length > intersections.second:
	containment = (
	GN_containment_types.throughgoing_bundle.value
	)

The tracks are checked if they intersect with the sphere, but there are sometimes tracks that decay after intersecting the sphere and before entering the hull. In this case, intersections.second will be negative, and we should not include this track in the extracted energy.

Additionally, sometimes none of the tracks that intersect the sphere make it to the detector. This can happen for single stopping muons from a CR. In this case, we run into even bigger problems after the for loop above, and one needs to check for that edge case and handle it the same way as if there were no tracks found at all.

[sevmag]

Getting the energies at certain track positions, like it is done here:

graphnet/src/graphnet/data/extractors/icecube/i3highesteparticleextractor.py

Line 322 in f3faff1

if MGtrack.get_energy(intersections.first) > EonEntrance:

can fail due to a MuonGun side assertion. We should catch these to avoid the conversion to crash.

Further appearances are:

graphnet/src/graphnet/data/extractors/icecube/i3highesteparticleextractor.py

Line 700 in f3faff1

track_energy = MGtrack.get_energy(intersections.first)

graphnet/src/graphnet/data/extractors/icecube/i3highesteparticleextractor.py

Line 333 in f3faff1

EonEntrance = MGtrack.get_energy(

[sevmag]

Looking at the following block in highest_energy_track

graphnet/src/graphnet/data/extractors/icecube/i3highesteparticleextractor.py

Lines 322 to 335 in f3faff1

	if MGtrack.get_energy(intersections.first) > EonEntrance:
	particle = track_particle
	closest_pos = np.array(
	[
	track.GetXc(),
	track.GetYc(),
	track.GetZc(),
	]
	)
	EonEntrance = MGtrack.get_energy(
	max(intersections.first, 0)
	)

Shouln't we already compute

MGtrack.get_energy(
    max(intersections.first, 0)
)

In the if statement? @Aske-Rosted