Questions about energy profile output and collision calculation

[Xiaozizhuo]

Hello dear developers, I am studying the code for this example
count_final_state_particles

Here I would like to raise two questions:
1.How to draw the energy distribution of each particle after 10,000 collisions on the basis of this code?

2.Why is there a count of the number of neutrinos of various flavors in the chart of the number of final particle states?

On the second question, my personal opinion is that the result of your defined code is not to show the total number of final particles produced after the collision, but the number of various mesons obtained after the collision plus the two total processes of decay, and the number of neutrinos to decay. Is my understanding correct?

[jncots]

1.How to draw the energy distribution of each particle after 10,000 collisions on the basis of this code?

Energy for particles of specific type can be obtained

pid = 2212 # proton or
pid = chromo.util.name2pdg("p")

for event in m(nevents):
     fs = event.final_state()
     protons_energies = fs.en[pid]

Using boost histogram:

# Energy histogram:
h_energy = bh.Histogram(
    bh.axis.Regular(50, 1e0, 1e3, transform=bh.axis.transform.log),
    bh.axis.IntCategory([], growth=True),
    bh.axis.Integer(0, len(models)),
)

mnames = []
for iModel, Model in enumerate(models):
    m = Model(kin, seed=1)
    mnames.append(m.name + " " + m.version)
    for event in m(nevents):
        fs = event.final_state()
        h.fill(fs.eta, fs.pid, iModel)
        h_energy.fill(fs.en, fs.pid, iModel) # <- add this

To plot use something like this:

from matplotlib import pyplot as plt
import numpy as np
from particle import Particle
from chromo.util import pdg2name

for im, mn in enumerate([mnames[0]]):
    print(mn)
    for iptype in range(len(h_energy.axes[1])):
        plt.stairs(h_energy[:,iptype, im].values(), h_energy.axes[0].edges, 
                label = pdg2name(h_energy.axes[1][iptype]))
        plt.xscale("log")
        plt.legend()

2.Why is there a count of the number of neutrinos of various flavors in the chart of the number of final particle states?

There is decay in the models. You can set stable/unstable particles individually like:

for pid in [111, 211, -211]:
    m.set_stable(pid)

or based on the decay time:

m.final_state_particles = chromo.util.select_long_lived(1e-20) # particles with lifetime > 1e-20 sec will be stable

[Xiaozizhuo]

Thank you for your help, now I have successfully implemented the operation of not letting particles decay. But the code that draws the energy distribution seems to be wrong.Here is the code I run.

import chromo
from matplotlib import pyplot as plt
import numpy as np
from particle import Particle
from chromo.util import pdg2name
from chromo.models import (
    DpmjetIII193,
    Phojet112,
    Sibyll23d,
)
from chromo.kinematics import CenterOfMass, GeV
import boost_histogram as bh
pid = 2212 # proton or
pid = chromo.util.name2pdg("p")
kin = CenterOfMass(1000 * GeV, "p", "p")
nevents = 10000
models = [Sibyll23d, Phojet112, DpmjetIII193]
h = bh.Histogram(
    bh.axis.Regular(20, -10, 10),
    bh.axis.IntCategory([], growth=True),
    bh.axis.Integer(0, len(models)),
)
h_energy = bh.Histogram(
    bh.axis.Regular(50, 1e0, 1e3, transform=bh.axis.transform.log),
    bh.axis.IntCategory([], growth=True),
    bh.axis.Integer(0, len(models)),
)
mnames = []
for iModel, Model in enumerate(models):
    m = Model(kin, seed=1)
    mnames.append(m.name + " " + m.version)
    for event in m(nevents):
        fs = event.final_state()
        h.fill(fs.eta, fs.pid, iModel)
        h_energy.fill(fs.en, fs.pid, iModel) # <- add this
for event in m(nevents):
     fs = event.final_state()
     protons_energies = fs.en[pid]
for im, mn in enumerate([mnames[0]]):
    print(mn)
    for iptype in range(len(h_energy.axes[1])):
        plt.stairs(h_energy[:,iptype, im].values(), h_energy.axes[0].edges, 
                label = pdg2name(h_energy.axes[1][iptype]))
        plt.xscale("log")
        plt.legend()
plt.show()

This is an error
IndexError: index 2212 is out of bounds for axis 0 with size 79

But when I change 2212 to the number that the individual particles represent, it still reports the same error

[jncots]

The part with pid was just an example how to get energy for specific particle type in simple case without boost_histogram. Otherwise

pid = 2212 # proton or
pid = chromo.util.name2pdg("p")
for event in m(nevents):
     fs = event.final_state()
     protons_energies = fs.en[pid]

is not needed

[Xiaozizhuo]

Thank you.It works!!!!