from chromo.models import Pythia8, EposLHC
from chromo.kinematics import FixedTarget, TeV, GeV
import numpy as np
from particle import literals as lp
from chromo.util import pdg2name
from rich.progress import trackWhen you look for a somewhat rare decay, you want to know how much luminosity is required to see N decays in your detector. For this one needs to compute the visual cross-section, this is the part of the production cross-section for the particle for which the decay products all fall into the acceptance of the detector.
We compute the visual cross-section here for the decay of Omega and Xi with chromo. We then go on and compute the expected rate of interactions for running LHCb in fixed-target mode.
Visual cross-section
We compute the inelastic cross-section for p-He collisions in fixed-target configuration where the proton has 13 TeV with chromo.
kin = FixedTarget(13*TeV, "p", "p")
# epos = EposLHC(kin)
gen = Pythia8(kin, seed=1)
cs = gen.cross_section()
*------------------------------------------------------------------------------------*
| |
| *------------------------------------------------------------------------------* |
| | | |
| | | |
| | PPP Y Y TTTTT H H III A Welcome to the Lund Monte Carlo! | |
| | P P Y Y T H H I A A This is PYTHIA version 8.308 | |
| | PPP Y T HHHHH I AAAAA Last date of change: 16 Nov 2022 | |
| | P Y T H H I A A | |
| | P Y T H H III A A Now is 31 Mar 2023 at 15:19:42 | |
| | | |
| | Program documentation and an archive of historic versions is found on: | |
| | | |
| | https://pythia.org/ | |
| | | |
| | PYTHIA is authored by a collaboration consisting of: | |
| | | |
| | Christian Bierlich, Nishita Desai, Leif Gellersen, Ilkka Helenius, Philip | |
| | Ilten, Leif Lonnblad, Stephen Mrenna, Stefan Prestel, Christian Preuss, | |
| | Torbjorn Sjostrand, Peter Skands, Marius Utheim and Rob Verheyen. | |
| | | |
| | The complete list of authors, including contact information and | |
| | affiliations, can be found on https://pythia.org/. | |
| | Problems or bugs should be reported on email at authors@pythia.org. | |
| | | |
| | The main program reference is C. Bierlich et al, | |
| | 'A comprehensive guide to the physics and usage of Pythia 8.3', | |
| | SciPost Phys. Codebases 8-r8.3 (2022) [arXiv:2203.11601 [hep-ph]] | |
| | | |
| | PYTHIA is released under the GNU General Public Licence version 2 or later.| |
| | Please respect the MCnet Guidelines for Event Generator Authors and Users. | |
| | | |
| | Disclaimer: this program comes without any guarantees. | |
| | Beware of errors and use common sense when interpreting results. | |
| | | |
| | Copyright (C) 2022 Torbjorn Sjostrand | |
| | | |
| | | |
| *------------------------------------------------------------------------------* |
| |
*------------------------------------------------------------------------------------*
cs.inelastic40.543016683445494Now we compute which fraction of inelastic events produce the following decays so that all charged final state decay products are in the LHCb acceptance, given by \(p > 2\) GeV/c and \(2 < \eta < 5\). \[ \Xi^- \to \Lambda \pi^- \; (99.887\,\%) \] \[ \Omega^- \to \Lambda K^- \; (67.8\,\%) \] We don’t need to take special care of the branching fractions. We get these reduction factors automatically if we select the specific decay we are interested in from the generator output.
apid_xi = abs(lp.Xi_minus.pdgid)
apid_omega = abs(lp.Omega_minus.pdgid)
apid_lambda = abs(lp.Lambda.pdgid)
apid_pi = abs(lp.pi_plus.pdgid)
apid_proton = abs(lp.proton.pdgid)
apid_K = abs(lp.K_plus.pdgid)
gen.set_stable(apid_xi, False)
gen.set_stable(-apid_xi, False)
gen.set_stable(apid_omega, False)
gen.set_stable(-apid_omega, False)
gen.set_stable(apid_lambda, False)
gen.set_stable(-apid_lambda, False)def is_visible_decay(ptot, eta, apid, children, par, apid1, apid2):
ch = children[par]
if len(ch) != 2:
return False
for i in ch:
if apid[i] != apid1 and apid[i] != apid2:
return False
if apid[i] == apid_lambda:
if not is_visible_decay(ptot, eta, apid, children, i, apid_proton, apid_pi):
return False
if ptot[i] < 2*GeV or eta[i] < 2 or eta[i] > 5:
return False
return True
n_inel = 100000
n_decays = {"Xi": 0, "Omega": 0}
for event in track(gen(n_inel), total=n_inel):
apid = np.abs(event.pid)
mask_xi = apid == apid_xi
mask_omega = apid == apid_omega
if np.any(mask_xi | mask_omega):
# collect children of each parent
children = {}
for i, par in enumerate(event.parents):
if par[0] == 0:
continue
children.setdefault(par[0] - 1, []).append(i)
ptot = event.p_tot
eta = event.eta
# investigate Xi decays
for par in np.arange(len(mask_xi))[mask_xi]:
if is_visible_decay(ptot, eta, apid, children, par, apid_lambda, apid_pi):
n_decays["Xi"] += 1
# investigate Omega decays
for par in np.arange(len(mask_omega))[mask_omega]:
if is_visible_decay(ptot, eta, apid, children, par, apid_lambda, apid_K):
n_decays["Omega"] += 1
print(n_decays, "out of", n_inel){'Xi': 122, 'Omega': 3} out of 100000sigma_vis = {k: v / n_inel * cs.inelastic for k, v in n_decays.items()}
print("visible cross-section")
for k, s in sigma_vis.items():
print(f"{k:10} {s:.4f} mb")visible cross-section
Xi 0.0495 mb
Omega 0.0012 mbLHCb in fixed-target mode
LHCb in fixed-target mode uses a gas target. We compute the event rate for the SMOG 1 configuration where the gas is injected directly into the LHCb Vertex Locator (VELO).
The average number of interactions per bunch crossing is \[ \mu = \sigma L N p \frac{N_A}{22400} \] where \(L = 60\) cm is the length of the detection region in the VELO, \(N=10^{11}\) is the proton bunch population in Run 2, \(p\) is gas pressure in bar, and \(N_A\) is Avogadro’s constant.
m = 1
b = 1e-28 * m ** 2
mb = 1e-3 * b
cm = 1e-2 * m
Pa = 1
bar = 1e5 * Pa
K = 1
kb = 1.380649e-23 # Boltzmann constant in J/K# cylinder volume in which a single proton would interact on average
cross_section = 200 * mb
detection_region_length = 60 * cm
volume = cross_section * detection_region_length
# number of collisions given by number of gas molecules per volume:
# pV = N k T -> N = p V / k T
p = 1e-10 * bar
T = 273 * K
n_collisions_per_proton = volume * p / (kb * T)
n_protons_per_bunch = 1e11
n_collisions = n_protons_per_bunch * n_collisions_per_proton
n_collisions0.0031837233037538114# Helium injection, 24/11/2022, target pressure on the VELO 1.2E-8mbar
p = 1.2e-8 * 1e-3