Question about how to obtain radius of innermost hit

[tmadlener]

Hi, I was looking into computing the radiusOfInnerMostHit because I need to use this variable for particle flow reconstruction. Since I have access to events generated with a previous software build that still has access to SiTracks_Refitted/SiTracks_Refitted.radiusOfInnermostHit, I was trying to confirm my piece of code below:

import math
import uproot

# load file
events = uproot.open("reco_p8_ee_tt_ecm365_100001.root")["events"]

# debug first event
iev = 0  
num_tracks = len(events["SiTracks_Refitted/SiTracks_Refitted.radiusOfInnermostHit"].array()[iev])

for itrack in range(num_tracks):
    
    val = events["SiTracks_Refitted/SiTracks_Refitted.radiusOfInnermostHit"].array()[iev][itrack]
    print(f"reading from branch; radiusOfInnermostHit of track {itrack} is: {val:.4f}")
    
    ############# now attempt to compute it using other branches
    # select the track states corresponding to itrack
    ibegin = events["SiTracks_Refitted/SiTracks_Refitted.trackStates_begin"].array()[iev][itrack]
    iend = events["SiTracks_Refitted/SiTracks_Refitted.trackStates_end"].array()[iev][itrack]

    num_track_states = iend - ibegin
    
    col = "_SiTracks_Refitted_trackStates"

    d0 = events[f"{col}/{col}.D0"].array()[iev][ibegin:iend]
    phi = events[f"{col}/{col}.phi"].array()[iev][ibegin:iend]
    refX = events[f"{col}/{col}.referencePoint.x"].array()[iev][ibegin:iend]
    refY = events[f"{col}/{col}.referencePoint.y"].array()[iev][ibegin:iend]    
    
    innermost_radius = float('inf')
    for ihit in range(num_track_states):
        x_closest = refX[ihit] + d0[ihit] * math.cos(phi[ihit])
        y_closest = refY[ihit] + d0[ihit] * math.sin(phi[ihit])
        
        radius = math.sqrt(x_closest**2 + y_closest**2)

        if radius < innermost_radius:
            innermost_radius = radius

    print(f"manually computed radiusOfInnermostHit of track {itrack} is: {innermost_radius:.4f}")
    break

which prints out:

reading from branch; radiusOfInnermostHit of track 0 is: 35.6223
manually computed radiusOfInnermostHit of track 0 is: 0.4091

Unfortunately, it seems that my piece of code is not correct because it does not match the value directly obtained from the branch SiTracks_Refitted/SiTracks_Refitted.radiusOfInnermostHit

Could you perhaps share how to build this variable from the other branches? Thanks!

Originally posted by @farakiko in #319 (comment)

[tmadlener]

We had quite a bit of discussion about this when we implemented the conversion from EDM4hep to LCIO for this particular piece of information in the corresponding PR: key4hep/k4EDM4hep2LcioConv#72 You can see some links in there to places where this is filled, and as you can also see, it's unfortunately not done consistently. Some places use a 2D distance, while others use a 3D distance. So depending on where this information is filled it is not easily possible to predict which distance was used in the end. Typically it is the 2D version IIRC.

Now in your code specifically, you do something like

x_closest = refX[ihit] + d0[ihit] * math.cos(phi[ihit])
y_closest = refY[ihit] + d0[ihit] * math.sin(phi[ihit])

whereas what we typically use is the x & y coordinate of the reference point directly, i.e.

x_closest = refX[ihit]
y_closest = refY[ihit]

[andresailer]

One should also only look at the track state atFirstHit, because the trackstate atIP will probably be smaller...
i.e.:
https://github.com/key4hep/k4EDM4hep2LcioConv/blob/019ba7030dbb7d3be71655a0b009dd3f8e1aa4fe/k4EDM4hep2LcioConv/src/k4EDM4hep2LcioConv.cpp#L171-L175

    if (state.location == edm4hep::TrackState::AtFirstHit) {
      const auto refP = state.referencePoint;
      return std::sqrt(refP.x * refP.x + refP.y * refP.y + use3D * refP.z * refP.z);
    }
  }

(The looping in the screenshot goes over track states, so calling the loop counter ihit is a bit confusing)

[farakiko]

Thanks for the feedback! If I do

x_closest = refX[ihit]
y_closest = refY[ihit]

I am still not recovering the values that I find in SiTracks_Refitted/SiTracks_Refitted.radiusOfInnermostHit

Is there something wrong in what I am looping over? Is there a mistake in the branch names?

[tmadlener]

You will have to first figure out which track state to use as simply looping over ihit and taking the first one is probably not giving you the one AtFirstHit (as Andre has pointed out above).

So you will also need to load the location branch for the track states and find the one (per track) where it matches edm4hep.TrackState.AtFirstHit.

[farakiko]

Gotcha! So I would pick the track state based on the location that coincides with edm4hep::TrackState::AtFirstHit and for that track state, compute the refX and refY to get the radius and that would be the radius of innermost hit?

If yes, where do I find edm4hep::TrackState::AtFirstHit? I can't seem to locate the branch in the rootfile

[andresailer]

It is the location branch of the trackstates.

edm4hep.TrackState.AtFirstHit is simply 2, if you don't want to import the whole edm4hep python bindings.

[tmadlener]

o I would pick the track state based on the location that coincides with edm4hep::TrackState::AtFirstHit and for that track state, compute the refX and refY to get the radius and that would be the radius of innermost hit?

Yes, exactly. But there is no computing involved any longer at that point, I think. refX and refY are the x & y coordinates that you want to use.

where do I find edm4hep::TrackState::AtFirstHit?

This is a constant value (effectively an enum). It's defined here

EDM4hep/edm4hep.yaml

Line 220 in fe5a540

static const int AtFirstHit = 2 ;\n

[farakiko]

Great! Thank you! Now the values match!!