Signal Muons Plot Issue

[A-A-Abdelhamid]

Here is the plot for signal muons Pts (status ==1 and is produced by a smuon vertex):

And the plot for final status muons Pts:

[A-A-Abdelhamid]

There is an issue with the plotting code, fix it using this:

import pyhepmc
import ROOT

def get_final_muon_descendant(particle):
    """
    Get the final state muon descendant of a given particle.

    Parameters
    ----------
    particle : pyhepmc.GenParticle
        The particle to check.
    """
    # Check if the particle itself is a final state muon
    if particle.status == 1 and abs(particle.pid) == 13:
        return particle
    # If the particle has an end vertex, check its descendants recursively
    elif particle.end_vertex:
        for p in particle.end_vertex.particles_out:
            final_muon = get_final_muon_descendant(p)
            if final_muon is not None:
                return final_muon
    # If the particle is not a final state muon and does not have an end vertex, return None
    return None

def fill_histogram_with_signal_muon_pt(hepmc_file):
    """
    Fill a histogram with the pt of signal muons in a HepMC file. 
    A signal muon is defined as a muon (13 or -13) that is produced by a decaying 2000013 or -2000013 
    and that is either a final state muon or decays into a final state muon.

    Parameters
    ----------
    hepmc_file : str
        The path to the HepMC file.
    """
    # Create a histogram
    hist = ROOT.TH1F("hist", "Muon transverse momentum;pt;Counts", 50, 0, 100)

    # Open the HepMC file
    with pyhepmc.open(hepmc_file) as f:
        # Loop over events in the file
        for event in f:
            # Loop over particles in the event
            for particle in event.particles:
                # Check if the particle is a muon produced by a decaying 2000013 or -2000013
                if abs(particle.pid) == 13 and particle.production_vertex and any(p.pid in [2000013, -2000013] for p in particle.production_vertex.particles_in):
                    # Get the final muon descendant of the muon
                    final_muon = get_final_muon_descendant(particle)
                    if final_muon is not None:
                        hist.Fill(final_muon.momentum.pt)

    # Draw the histogram
    canvas = ROOT.TCanvas("canvas")
    hist.Draw()
    canvas.SaveAs("histogram.png")

# Use the function
fill_histogram_with_signal_muon_pt('tag_1_pythia8_events.hepmc')

[A-A-Abdelhamid]

The current code works as intended:

def get_final_muon_descendant(particle):

    #Get the final state muon descendant of a given particle.
    # Check if the particle itself is a final state muon
    if particle.status == 1 and abs(particle.pid) == 13:
        return particle
    # If the particle has an end vertex, check its descendants recursively
    elif particle.end_vertex:
        for p in particle.end_vertex.particles_out:
            final_muon = get_final_muon_descendant(p)
            if final_muon is not None:
                return final_muon
    # If the particle is not a final state muon and does not have an end vertex, return None
    return None

Then call it inside a for particle loop

for event in f:
      for particle in event.particles:
        momentum =particle.momentum
        pt = momentum.pt()

       # Check if the particle is a muon produced by a decaying 2000013 or -2000013
        if abs(particle.pid) == 13 and particle.production_vertex and any(p.pid in [2000013, -2000013] for p in particle.production_vertex.particles_in):
            # Get the final muon descendant of the muon
            final_muon = get_final_muon_descendant(particle)
            if final_muon is not None:
                histSignalPt.Fill(final_muon.momentum.pt())

Now we get the correct signal muons pt plot:

This does not impact the final status muon pt plot:

[A-A-Abdelhamid]

Looking at one of the events an (anti) muon with status ==1 but not produced by a smuon decay:

It is produced by a meson decay, its pt is ~ 1 GeV, total energy of ~12 GeV