analysis_toolkit: Candidate Selection and Event Inspection

This module lives under the experimental package. Its API may change significantly between FairShip releases. It is a work-in-progress and should not be used in production analyses without review.

The experimental.analysis_toolkit module provides two high-level classes for reconstructed-event analysis in SHiP. selection_check encapsulates the full standard pre-selection pipeline: it loads the detector geometry from a ROOT file, applies impact parameter, fiducial volume, track quality, and timing cuts on HNL decay candidates, and can render a colour-coded table of per-cut pass/fail results. event_inspector provides a complementary view at the Monte Carlo truth level, printing a formatted table of all MCTrack entries per event for rapid visual debugging.

import ROOT
from experimental import analysis_toolkit

# Open files
f        = ROOT.TFile.Open("ship.conical.Pythia8-TGeant4.root", "read")
tree     = f["cbmsim"]
tree.AddFriend("ship_reco_sim", "ship.conical.Pythia8-TGeant4_rec.root")
geo_file = ROOT.TFile.Open("geofile_full.conical.Pythia8-TGeant4.root", "read")

# Initialise toolkit objects
selection = analysis_toolkit.selection_check(geo_file)
inspector = analysis_toolkit.event_inspector()

for event_nr, event in enumerate(tree):
    selection.access_event(event)          # bind current event tree
    inspector.dump_event(event, mom_threshold=0.5)

    for candidate in event.Particles:
        passed = selection.preselection_cut(candidate, IP_cut=250, show_table=True)
        if passed:
            print(f"Event {event_nr}: candidate passes pre-selection")

`selection_check`

selection_check wraps all geometry-aware per-candidate calculations and pre-selection cuts. It must be initialised with a ROOT geometry file so it can set up the TGeoManager and read the ShipGeo configuration.

Constructor

selection_check(geo_file)

geo_file

ROOT.TFile

required

An open ROOT.TFile containing both the FAIRGeom geometry manager and the ShipGeo configuration object (in either JSON or pickle format). The geometry manager is retrieved via geo_file.Get("FAIRGeom") and the configuration via load_from_root_file.On construction the veto geometry YAML is also loaded from $FAIRSHIP/geometry/:

veto_config_helium.yaml when ship_geo.DecayVolumeMedium == "helium"
veto_config_vacuums.yaml when ship_geo.DecayVolumeMedium == "vacuums"

Attributes set by __init__:

Attribute	Type	Description
`geometry_manager`	`TGeoManager`	ROOT geometry manager loaded from the geo file
`ship_geo`	`Config`	Full geometry configuration object
`veto_geo`	`AttrDict`	Veto detector geometry parameters from YAML

`access_event`

access_event(tree) -> None

Bind the current event tree to the selection_check instance. This must be called once per event (or once per tree if iterating with for event in tree) before invoking any per-candidate method, since those methods read branches such as FitTracks, fitTrack2MC, and strawtubesPoint from self.tree.

tree

ROOT.TTree | event

required

The current event object obtained from iterating the cbmsim tree (with the reconstruction friend attached).

`define_candidate_time`

define_candidate_time(candidate) -> float

Estimate the time of the decay vertex for a candidate by back-propagating the timing of straw-tube MC hits from the first two tracking stations. For each straw hit associated with a daughter track, the propagation time from the decay vertex to the hit position is subtracted (assuming straight-line travel at the particle’s velocity), then the event start time t0 from the ShipEventHeader is added.

candidate

FairParticle / reconstructed particle

required

The candidate decay vertex object. Must expose GetVertex(), GetDaughter(i), and Momentum() methods, as provided by the SHiP reconstruction output.

return

float

Estimated time at the decay vertex in nanoseconds (FairShip internal time units), averaged over all contributing straw hits from daughter tracks.

`impact_parameter`

impact_parameter(candidate) -> float

Compute the distance of closest approach of the candidate’s momentum vector to the target position (0, 0, ship_geo.target.z0) — the standard impact parameter used to reject combinatorial background.

candidate

reconstructed particle

required

The candidate object. Vertex and four-momentum are read via GetVertex() and Momentum().

return

float

Impact parameter in cm.

`dist_to_innerwall`

dist_to_innerwall(candidate) -> float | int

Find the minimum distance in the XY plane from the candidate decay vertex to the inner wall of the decay vessel, by ray-casting in 8 equally-spaced azimuthal directions and querying TGeoManager::FindNextBoundary. Returns 0 if the vertex is outside the decay volume or if the wall distance exceeds 100 m.

candidate

reconstructed particle

required

Candidate whose vertex is tested.

return

float | int

Minimum distance to the inner wall in cm, or 0 if outside the volume or wall distance exceeds 100 m.

`dist_to_vesselentrance`

dist_to_vesselentrance(candidate) -> float

Return the z-distance of the candidate decay vertex from the entrance plane of the decay vessel, defined as vertex.Z() - veto_geo.z0.

candidate

reconstructed particle

required

Candidate whose vertex z-coordinate is tested.

return

float

Distance in cm from the vessel entrance. Positive values are downstream of the entrance.

`is_in_fiducial`

is_in_fiducial(candidate) -> bool

Check whether the candidate decay vertex lies within the fiducial decay volume. The vertex must be:

Downstream of the veto entrance (vertex.Z() > veto_geo.z0)
Upstream of Track Station 1 (vertex.Z() < ship_geo.TrackStation1.z)
Inside a volume whose ROOT geometry name starts with "DecayVacuum_"

candidate

reconstructed particle

required

Candidate to test.

return

bool

True if the decay vertex is within the defined fiducial volume.

`nDOF`

nDOF(candidate) -> np.ndarray

Return the number of degrees of freedom from the track fits of the two daughter tracks.

candidate

reconstructed particle

required

Candidate providing daughter track indices via GetDaughter(0) and GetDaughter(1).

return

numpy.ndarray

Array of shape (2,) containing the rounded nDOF for the first and second daughter tracks, read from FitTracks[i].getFitStatus().getNdf().

`chi2nDOF`

chi2nDOF(candidate) -> np.ndarray

Return the reduced χ²/nDOF from the track fits of the two daughter tracks.

candidate

reconstructed particle

required

Candidate providing daughter track indices.

return

numpy.ndarray

Array of shape (2,) with chi2 / ndf for each daughter track.

`daughtermomentum`

daughtermomentum(candidate) -> np.ndarray

Return the fitted momentum magnitude of each daughter track.

candidate

reconstructed particle

required

Candidate providing daughter track indices.

return

numpy.ndarray

Array of shape (2,) with momentum magnitude in GeV for each daughter, read from FitTracks[i].getFittedState().getMom().Mag().

`invariant_mass`

invariant_mass(candidate) -> float

Return the reconstructed invariant mass of the candidate.

candidate

reconstructed particle

required

Candidate exposing GetMass().

return

float

Invariant mass in GeV.

`DOCA`

DOCA(candidate) -> float

Return the Distance of Closest Approach between the two daughter tracks at the decay vertex.

candidate

reconstructed particle

required

Candidate exposing GetDoca().

return

float

DOCA in cm.

`preselection_cut`

preselection_cut(candidate, IP_cut: float = 250, show_table: bool = False) -> bool

Apply the full standard pre-selection to a candidate and return True if all cuts pass. This is the umbrella method that calls all individual selection methods. The applied cuts are:

Criterion	Cut
Number of candidates in event	`== 1`
Fiducial volume	`is_in_fiducial == True`
Distance to inner wall	`> 5 cm`
Distance to vessel entrance	`> 100 cm`
Impact parameter	`< IP_cut` (in cm)
DOCA	`< 1 cm`
nDOF (both daughters)	`> 25`
χ²/nDOF (both daughters)	`< 5`
Daughter momentum (both)	`> 1 GeV`

candidate

reconstructed particle

required

The candidate to evaluate.

IP_cut

float

default:"250"

Maximum allowed impact parameter in cm. Default is 250 cm.

show_table

bool

default:"False"

When True, prints a formatted tabulate grid table to stdout showing each criterion’s measured value, cut threshold, and pass/fail status (ANSI-coloured green/red).

return

bool

True if the candidate passes all pre-selection cuts, False otherwise.

for candidate in event.Particles:
    # Quick pass/fail
    if selection.preselection_cut(candidate):
        print("passed")

    # With detailed per-cut printout
    selection.preselection_cut(candidate, IP_cut=200, show_table=True)

`event_inspector`

event_inspector provides a Monte Carlo truth viewer. It initialises the PDG database and registers the HNL particle type so that HNL tracks display with a meaningful name.

Constructor

event_inspector()

Initialises ROOT.TDatabasePDG.Instance() and calls pythia8_conf.addHNLtoROOT() to register HNL PDG entries.

`dump_event`

dump_event(event, mom_threshold: float = 0) -> None

Print a formatted table of all MCTrack entries in the event using tabulate. Each row shows the track index, particle name (from PDG database, or "----" if unknown), PDG code, mother track ID, momentum components (Px, Py, Pz) in GeV/c, start vertex (x, y, z) in metres, Geant4 process name, and event weight.

event

ROOT.TTree event

required

The current event object, expected to contain a MCTrack branch.

mom_threshold

float

default:"0"

Minimum momentum in GeV for a track to be included in the table. Tracks with P < mom_threshold are skipped. Useful for suppressing low-momentum secondaries.

inspector = analysis_toolkit.event_inspector()

for event_nr, event in enumerate(tree):
    print(f"\n=== Event {event_nr} ===")
    inspector.dump_event(event, mom_threshold=0.5)
    # Prints a table with columns:
    # # | particle | pdgcode | mother_id | Momentum [Px,Py,Pz] (GeV/c) |
    #   | StartVertex[x,y,z] (m) | Process | GetWeight()

Full example from `examples/analysis_example.py`

The following is the complete analysis workflow taken directly from examples/analysis_example.py:

#!/usr/bin/env python
"""Example script for usage of the analysis_toolkit for signal selection."""

from argparse import ArgumentParser
import ROOT
import rootUtils as ut
from experimental import analysis_toolkit


def main():
    parser = ArgumentParser(description=__doc__)
    parser.add_argument("-f", "--simfile",  default="ship.conical.Pythia8-TGeant4.root")
    parser.add_argument("-r", "--recofile", default="ship.conical.Pythia8-TGeant4_rec.root")
    parser.add_argument("-g", "--geofile",  default="geofile_full.conical.Pythia8-TGeant4.root")
    options = parser.parse_args()

    # Open simulation and reconstruction files
    f    = ROOT.TFile.Open(options.simfile, "read")
    tree = f["cbmsim"]
    tree.AddFriend("ship_reco_sim", options.recofile)

    geo_file  = ROOT.TFile.Open(options.geofile, "read")
    selection = analysis_toolkit.selection_check(geo_file)
    inspector = analysis_toolkit.event_inspector()

    # Book histograms for pre-selection observables
    hist_dict = {}
    ut.bookHist(hist_dict, "event_weight", "Event weight", 100, 100, 100)
    ut.bookHist(hist_dict, "candidate_time", "candidate time @ decay vertex; ns", 200, 100, 100)
    ut.bookHist(hist_dict, "impact_parameter", "Impact parameter; cm", 200, 100, 100)
    ut.bookHist(hist_dict, "dist_to_innerwall", "Distance to inner wall; cm", 200, 100, 100)
    ut.bookHist(hist_dict, "dist_to_vesselentrance", "Distance to Decay Vessel Entrance; cm", 200, 100, 100)
    ut.bookHist(hist_dict, "inv_mass", "Invariant mass; GeV", 200, 100, 100)
    ut.bookHist(hist_dict, "DOCA", "Distance of closest approach; cm", 200, 100, 100)
    ut.bookHist(hist_dict, "len_Particles", "len(tree.Particles); Number of candidates per event", 200, 100, 100)
    ut.bookHist(hist_dict, "d_mom", "momentum of daughters; d1 (GeV); d2 (GeV)", 200, 100, 100, 200, 100, 100)
    ut.bookHist(hist_dict, "nDOF", "nDOF; d1; d2", 200, 100, 100, 200, 100, 100)
    ut.bookHist(hist_dict, "chi2nDOF", "chi2nDOF; d1; d2", 200, 100, 100, 200, 100, 100)

    for event_nr, event in enumerate(tree):
        selection.access_event(event)
        if len(event.Particles) == 0:
            continue

        print(f"Event{event_nr}:")
        inspector.dump_event(event, mom_threshold=0.5)

        event_weight = event.MCTrack[2].GetWeight()
        hist_dict["event_weight"].Fill(event_weight)

        for candidate_id_in_event, signal in enumerate(event.Particles):
            hist_dict["candidate_time"].Fill(selection.define_candidate_time(signal))
            hist_dict["impact_parameter"].Fill(selection.impact_parameter(signal))
            hist_dict["dist_to_innerwall"].Fill(selection.dist_to_innerwall(signal))
            hist_dict["dist_to_vesselentrance"].Fill(selection.dist_to_vesselentrance(signal))
            hist_dict["DOCA"].Fill(selection.DOCA(signal))
            hist_dict["inv_mass"].Fill(selection.invariant_mass(signal))
            hist_dict["len_Particles"].Fill(len(tree.Particles))
            hist_dict["d_mom"].Fill(*selection.daughtermomentum(signal))
            hist_dict["nDOF"].Fill(*selection.nDOF(signal))
            hist_dict["chi2nDOF"].Fill(*selection.chi2nDOF(signal))

            IP_cut = 250
            preselection_flag = selection.preselection_cut(signal, IP_cut=IP_cut, show_table=False)

            if preselection_flag:
                print(f"Event:{event_nr} Candidate_index: {candidate_id_in_event} <--passes the pre-selection\n\n")

    ut.writeHists(hist_dict, "preselectionparameters.root")


if __name__ == "__main__":
    main()

Python Modules

Configuration

analysis_toolkit: Candidate Selection and Event Inspection

`selection_check`

Constructor

`access_event`

`define_candidate_time`

`impact_parameter`

`dist_to_innerwall`

`dist_to_vesselentrance`

`is_in_fiducial`

`nDOF`

`chi2nDOF`

`daughtermomentum`

`invariant_mass`

`DOCA`

`preselection_cut`

`event_inspector`

Constructor

`dump_event`

Full example from `examples/analysis_example.py`

Build docs developers (and LLMs) love

Python Modules

Configuration

Documentation Index

​selection_check

​Constructor

​access_event

​define_candidate_time

​impact_parameter

​dist_to_innerwall

​dist_to_vesselentrance

​is_in_fiducial

​nDOF

​chi2nDOF

​daughtermomentum

​invariant_mass

​DOCA

​preselection_cut

​event_inspector

​Constructor

​dump_event

​Full example from examples/analysis_example.py

Build docs developers (and LLMs) love

`selection_check`

Constructor

`access_event`

`define_candidate_time`

`impact_parameter`

`dist_to_innerwall`

`dist_to_vesselentrance`

`is_in_fiducial`

`nDOF`

`chi2nDOF`

`daughtermomentum`

`invariant_mass`

`DOCA`

`preselection_cut`

`event_inspector`

Constructor

`dump_event`

Full example from `examples/analysis_example.py`