Complete decay chain issue

[A-A-Abdelhamid]

Using this code

import pyhepmc
import ROOT
from collections import Counter

def get_decay_products(particle, event_num):
    """
    Get the decay products of a given particle.
    """
    proper_decay_products = []
    other_particles = []

    # Check if the particle itself is a final state muon or neutrino
    if particle.status == 1 and abs(particle.pid) in [13, 12, -13, -12]:
        proper_decay_products.append((particle.pid, particle.id))
    # If the particle has an end vertex, check its descendants recursively
    elif particle.end_vertex:
        for p in particle.end_vertex.particles_out:
            if abs(p.pid) in [13, 12, -13, -12]:
                proper, other = get_decay_products(p, event_num)
                proper_decay_products.extend(proper)
                other_particles.extend(other)
            elif p.pid!=22:
                other_particles.append((p.pid, p.id, event_num))
    return proper_decay_products, other_particles

def track_decay_chain(hepmc_file):
    """
    Track the decay chain of muons and neutrinos that originate from a particle with abs pid 2000013.
    """
    proper_decay_products = []
    other_particles = []

    # Open the HepMC file
    with pyhepmc.open(hepmc_file) as f:
        # Loop over events in the file
        for i, event in enumerate(f):
            # Loop over particles in the event
            for particle in event.particles:
                # Check if the particle is a muon or neutrino produced by a decaying 2000013 or -2000013
                if abs(particle.pid) in [13, 12, -13, -12] and particle.production_vertex and any(p.pid in [2000013, -2000013] for p in particle.production_vertex.particles_in):
                    # Get the decay products of the particle
                    proper, other = get_decay_products(particle, i)
                    proper_decay_products.extend(proper)
                    other_particles.extend(other)

    # Count the frequencies of each particle PID
    proper_decay_counter = Counter([pid for pid, _ in proper_decay_products])
    other_particles_counter = Counter([pid for pid, _, _ in other_particles])

    print("Proper decay products:", proper_decay_counter)
    print("Other particles:", other_particles_counter)

track_decay_chain('tag_1_pythia8_events.hepmc')

I got the output:

Proper decay products: Counter({13: 20000, 12: 20000, -13: 20000, -12: 20000})
Other particles: Counter()

So, by excluding the photons from the chain, everything works as expected.

Now by changing elif p.pid!=22: to anelse:

Proper decay products: Counter({13: 20000, 12: 20000, -13: 20000, -12: 20000})
Other particles: Counter({22: 24498})

Thus it is clear that the other unexpected particles come from photons, since the the code only check muons and neutrinos recursively, but stop at any other particle after appending it to "other particles"

[A-A-Abdelhamid]

Now checking the entire decay chain recursively for every particle (not only muons and neutrinos) by changing the code for "get_decay_products" function to:

def get_decay_products(particle, event_num):
    """
    Get the decay products of a given particle.

    """
    proper_decay_products = []
    other_particles = []

    # Check if the particle itself is a final state muon or neutrino
    if particle.status == 1 and abs(particle.pid) in [13, 12, -13, -12]:
        proper_decay_products.append((particle.pid, particle.id))
    elif particle.status == 1:
        other_particles.append((particle.pid, particle.id, event_num))
    # If the particle has an end vertex, check its descendants recursively
    if particle.end_vertex:
        for p in particle.end_vertex.particles_out:
            proper, other = get_decay_products(p, event_num)
            proper_decay_products.extend(proper)
            other_particles.extend(other)
    return proper_decay_products, other_particles

We get this output:

Proper decay products: Counter({13: 134849, -13: 134840, -12: 134688, 12: 134679})
Other particles: Counter({22: 190298, -211: 956, 211: 929, 11: 535, -11: 533, 321: 132, 
2212: 111, -2212: 111, 
130: 92, -321: 87, -14: 19, -2112: 10, 2112: 10, -16: 5, 16: 5, 14: 3})

Please notice that this code overcounts particles, I'm looking into it, but I looked at one of the desired events. I already verified that such particles only appear when children particles of photons are considered.

[A-A-Abdelhamid]

Currently looking at the events of those unexpected particles

[A-A-Abdelhamid]

# Find the first non photon PID in other_particles, which is the first particle that is out of a photon end_vertix
    first_non_22 = next((pid, id, event_num) for pid, id, event_num in other_particles if pid != 22)

Output:

First non-22 PID in 'other_particles': PID = 11, ID = 1399, Event number = 11

[A-A-Abdelhamid]

Event 11

[A-A-Abdelhamid]

Here is what is happening:

The photon created a pair of e- e+ with different energies.

[A-A-Abdelhamid]

This is the code used to print the event diagram with the particle of interest and its vertex position

import pyhepmc
from pyhepmc.view import savefig

filename = "tag_1_pythia8_events.hepmc"
# pyhepmc.open can read most HepMC formats using auto-detection
with pyhepmc.open(filename) as f:
    # Loop over events in the file
    for i, event in enumerate(f):
        # Save a plot of the event
        if event.event_number == 11:
          savefig(event, f"event{i}.svg")
          savefig(event, f"event{i}.png")
          savefig(event, f"event{i}.pdf")
          for particle in event.particles:
            if particle.id ==1399 and particle.pid==11:
              ver=particle.production_vertex
              print(ver)

However, it creates diagrams for every single event with "11" in its event number, but that is not a concern at the moment.

[A-A-Abdelhamid]

I checked one of the events where we get unexpected particles, this particular one was an electron coming from an e- e+ pair production by a photon. All unexpected particles in the chain come from a vertex with a photon coming in. The process is detailed above for your reference. Thanks! @trholmes

[trholmes]

I don't totally understand from looking at your diagram. Yes the one of your muons radiates two photons, and the photons create e+e- pairs, but I don't see any other particles coming from those chains. Can you clarify?

[A-A-Abdelhamid]

This is an event where unexpected particles appear, the e+ e- themselves are not expected particles from our decay chain (smuon ---> mu + neutrino). It just happened that the photon in this event was creating two final status electrons, but if I pick another event, I could for example see a photon creating pi+ pi- I hope this clarifies what I meant.

[A-A-Abdelhamid]

e-e+ was not the best example. Here is another example, the most observed unexpected particle in our chain is -211 (pi-)

Here is a pi- in a photon chain (where the photon itself is from a muon coming from smuon) PID = -211, ID = 552, Event number = 868 production vertex: (FourVector(-10.8, -21.9, 22.9, 46.4)) Event #868

Here is the decay chain:

The same picture but clearer (the chain of the photon with 58 GeV)

This one has a lot of unexpected particles compared to photon -----> e- e+ e-e+ are also not expected from our smuon decay chain and I should have mentioned that.

I hope this clarifies everything.