SoftKiller should return PseudoJets with the correct `cluster_hist_index`

[graeme-a-stewart]

Currently SoftKiller is copying out PseudoJet objects verbatim using filter:

JetReconstruction.jl/src/SoftKiller.jl

Line 176 in 1270f05

reduced_event = filter(jet -> pt2(jet) >= pt2cut, event)

There are two issues here:

1. This is quite expensive, but at the moment I don't see a way to get around this: we need a real vector of clusters to then initiate jet reconstruction, which will be the next step in the chain (if we were to return, e.g., references, then we'd need to construct the clusters vector at a later stage anyway), so we're just paying a necessary cost up-front that needs to be paid.
2. The cluster_hist_index values are wrong, because they come from the cluster order as read from read_final_state_particles.

The consequence of 2 is that the user needs to recreate the reduced_event again with the correct sequential values for the cluster_hist_index. e.g. as done here:

JetReconstruction.jl/examples/softkiller/softkiller_plots.jl

Lines 88 to 92 in 1270f05

	new_pseudo_jet = PseudoJet(JetReconstruction.px(pseudo_jet),
	JetReconstruction.py(pseudo_jet),
	JetReconstruction.pz(pseudo_jet),
	JetReconstruction.energy(pseudo_jet);
	cluster_hist_index = i)

).

This also means that the example given in the documentation is wrong as the reconstruction will crash when the @assert which checks for sequential cluster_hist_index will fail.

So we're paying the price of two copies, when we only really need one.

So, I propose that instead of using filter we loop over the event's clusters and create new clusters, with the correct linear ordering for cluster_hist_index.

N.B. for reasons of efficiency it is important to use the fully specified constructor for PseudoJet, e.g.,

new_cluster = PseudoJet(c.px, c.py, cpz, c.E, the_new_cluster_hist_index, c. _pt2, c._inv_pt2, c._rap, c._phi)

This will avoid the expensive recalculation of the internally cached values.

[emadmtr]

Hi @graeme-a-stewart I cheged the code from using the filter to creating a new PseudoJet in a for loop, should I open a new PR to commit the changes? Also I have a question, in the original FastJet implementation they use the class Selector to remove any particles with rapidity smaller then -5 and greated than 5, there I beleive I also use a filter (in method select_ABS_RAP_max) do you thin k I should also change it, it shouldn't be affecting the SoftKiller efficiency itself but it might be overall a bad choce of approaching the problem.

[emadmtr]

Hi @graeme-a-stewart, so sorry for bothering you, I know you are probably very busy, but I just wanted to follow up on my last comment.

[graeme-a-stewart]

Hello @emadmtr - sorry for the late reply, we're in the middle of a relocation, so many non-work things happening!

Yes, please open a PR with the changes. Were you able to check if this works just as fast as the filter() or if we're playing a small price for this?

Re. the rapidity filter, in general it's good to reproduce the Fastjet behaviour. However, it would seem useful to make this configurable, with a default range of ±5.

Feel free to poke me again if there's something that you're waiting on!

Graeme