Is there a way to get the constituent indices of subjets with array-based interface?

[jmduarte]

Hi all, excellent work on this python interface!

I was wondering if there's a way to get the constituent indices of the subjets using the array-based interface. For example, I have an awkward array of constituents corresponding to a jet clustered with anti-kT R=0.8. I can recluster them with the same algorithm (resulting in 1 jet as expected), and then get subjets with exclusive_subjets(). But I'm not sure if there's a way to get the constituent indices that correspond to those subjets.

It seems there may be a way to do it with the object-based approach (because each PseudoJet object can return its own constituents), but then it has to be done one jet/event at a time. Is this correct?

Thank you!

[jpivarski]

If I remember right, constituents was one of the member functions of the vectorized Coasters object (the same object you get exclusive subjets from. They're direct links to them, so if the input particles had non-vector attributes (like "charge" or "isolation"), then they're carried through to the constituents.

[jpivarski]

(Now that I'm on a computer...) The actual constituents are in a method called ClusterSequence.constituents, but that should be the same as array-slicing the input data with constituent_index.

In both cases, these functions are supposed to have the same length as the inclusive_jets, so this doesn't look right. That's why the output is a doubly nested list. (Just confirmed: page 11 of @aryan26roy's talk here.)

I suspect you'll find something odd in the cluster.constituent_index().layout. Do you find the same thing with antikt_algorithm? Most of the tests and demos were done with that.

So, this project needs a maintainer, a point that I raised in IRIS-HEP Slack #pyjet-and-beyond and elsewhere. (I'll invite you there.) I wanted to try testing it, but my Linux computer broke and I'm using a Mac, it has an M2 chip, and a wheel hasn't been compiled for that case, the compilers I've installed fail on first try... As it is, a few physics groups do find this interface useful, but the package in its current state has some rough edges. The original developer (@aryan26roy) is not a physicist—this kind of "last 10%" clean-up needs to come from physicist users who know the use-cases.

[jmduarte]

Yes, it's the same with anti-kT.

To clarify, the length of cluster.constituent_index() is the same as cluster.inclusive_jets()! They both are length 1 in the innermost dimension because there is exactly 1 jet. This is probably as it should be because I'm only taking the constituents of an existing AK8 jet and re-clustering it.

The layout is as follows:

<ListOffsetArray64>
    <offsets><Index64 i="[0 1]" offset="0" length="2" at="0x562637244da0"/></offsets>
    <content><ListOffsetArray64>
        <offsets><Index64 i="[0 55]" offset="0" length="2" at="0x56263721b5c0"/></offsets>
        <content><NumpyArray format="i" shape="55" data="0 1 2 3 4 ... 50 51 52 53 54" at="0x562636b89630"/></content>
    </ListOffsetArray64></content>
</ListOffsetArray64>

So I think everything is working correctly. My problem is I don't know if it's possible with the current interface to find out which constituents belong to the 3 exclusive subjets of the 1 reconstructed jet.

For testing, I updated my gist to be a notebook that you can run in Colab!
https://gist.github.com/jmduarte/e9f981a6a22834ba7b3d9a675a1b915a

Re: maintenance, I completely understand. If my students/postdocs/myself can help to fix some of these issues "last 10%" issues with PRs, we certainly are happy to! Again, I'm not even sure if this is an "issue" as much as a "feature request"

[jpivarski]

Ah, I see. I never quite understood the difference between the "inclusive" and "exclusive" options. I think of jet-finding as hierarchical/agglomerative clustering. For the exclusive_subjets, you specify a number of clusters (nsub=3); is that k-means? And is that subdividing the already-found cluster or is it re-clustering the original particles?¹

Thanks for the Google Colab version. It didn't let me import; I had to copy-paste the code blocks (which is something I could have done from the beginning, had I thought of Colab—I was thinking about launching AWS or Binder, which is more set-up), but it worked. I was able to verify that the number of constituents follows the number of inclusive_jets by varying the ΔR:

jetdef = fastjet.JetDefinition(fastjet.cambridge_algorithm, 0.3)
cluster = fastjet.ClusterSequence(array, jetdef)
jets = cluster.inclusive_jets()
ak.num(jets), ak.num(cluster.constituents())

(<Array [5] type='1 * int64'>, <Array [5] type='1 * int64'>)

There are a number of things that need to be done (https://github.com/scikit-hep/fastjet/issues) and I know of two user groups: yours and the SUEP group at MIT (Chad Freer and @chrispap95). Maybe we should get together on a call, assess what needs to be done, find an equitable way to share the burden, and how I can help guide the process/answer questions about the code?

Footnotes

1. This may be illustrating the point I made earlier. @aryan26roy is a computer scientist and although I'm a physicist, most of my actually-doing-physics-work was at a 10 GeV e⁺e⁻ collider (no jets) and I did a muon analysis at the LHC. Neither of us know a lot about jet-finding; we were just passing through the FastJet (capital F, J) interface in an array-oriented way. The question here is whether FastJet has a "constituents" function for exclusive_subjets. ↩

[jmduarte]

@jpivarski sure, I'd be happy to have a call and see if I can enlist some students to help

[jpivarski]

(We followed up by email.)

[aryan26roy]

Hi @jmduarte, sorry for the late reply. You are right, cluster.constituent_index() returns the constituents of the inclusive jets. As for your problem, in the C++ version, I presume that once you get the vector of PseudoJet objects from the cluster.exclusive_subjets() function call, you would go through each of them doing the PseudoJet.constituents() function call on each object to get what you want?

If so, there is no way of doing that with this interface right now. However, it can be easily implemented the same way cluster.constituent_index() was implemented for inclusive jets.

As @jpivarski mentioned, we were not too familiar with jet-finding, hence we only implemented the already available functions (and a few derivative ones) to get the first working version. These few implementations can easily serve as templates to add whatever functionality a user wants without much trouble.

[jmduarte]

Thanks @aryan26roy, this is what I figured. I just wasn't sure if there was some "trick" that can do it.

For anyone else following, here's the solution I landed on, based on looping over the pseudojets constituents and checking for equality with the original constituents: https://gist.github.com/jmduarte/e9f981a6a22834ba7b3d9a675a1b915a