@@ -155,30 +155,40 @@ def total_track_energy(
155155 track_ids = np .array ([track .id for track in track_list ])
156156
157157 total_tracks = len (track_list )
158+ exclude_ids = set ()
158159 while len (track_list ) > 0 :
159160 track = track_list [0 ]
160161
161162 # Find distance to entrance and exit from sampling volume
162163 intersections = self .hull .surface .intersection (
163164 track .pos , track .dir
164165 )
165- # Get the corresponding energies
166166
167+ # Check if the track actually enters the volume
168+ if not (
169+ np .isfinite (intersections .first )
170+ and (
171+ intersections .first
172+ < frame [self .mctree ].get_particle (track .id ).length
173+ )
174+ ):
175+ track_list = track_list [1 :]
176+ track_ids = track_ids [1 :]
177+ frame [self .mctree ].erase (track .id )
178+ self ._logger .debug (
179+ f"Remaining tracks: { len (track_list )} /{ total_tracks } "
180+ )
181+ continue
182+
183+ # Get the corresponding energies
167184 e0 = track .get_energy (intersections .first )
168185 e1 = track .get_energy (intersections .second )
169186
170187 # Accumulate
171188 e_deposited += e0 - e1
172189 e_entrance += e0
173- if self .daughters :
174- assert e_entrance <= sum (
175- [p .energy for p in primaries ]
176- ), "Energy on entrance is greater than primary energy"
177- assert e_deposited <= sum (
178- [p .energy for p in primaries ]
179- ), "Energy deposited is greater than primary energy"
180190 # erase particle and children
181- exclude_ids = set (
191+ exclude_ids . update (
182192 [c .id for c in frame [self .mctree ].children (track .id )]
183193 )
184194 exclude_ids .add (track .id )
@@ -190,6 +200,14 @@ def total_track_energy(
190200 self ._logger .debug (
191201 f"Remaining tracks: { len (track_list )} /{ total_tracks } "
192202 )
203+ # Sanity check ensuring no double counting
204+ if self .daughters :
205+ assert e_entrance <= sum (
206+ [p .energy for p in primaries ]
207+ ), "Energy on entrance is greater than primary energy"
208+ assert e_deposited <= sum (
209+ [p .energy for p in primaries ]
210+ ), "Energy deposited is greater than primary energy"
193211 return e_entrance , e_deposited
194212
195213 def total_cascade_energy (
0 commit comments