Setting up simulation of p-A collisions

[michael-pitt]

Hello,

In p-A collisions, hadrons have asymmetric energy (E_A = Z*E_P != E_P). If kinematics are defined using the CenterOfMass class, then the results must to be boosted to the lab frame. Is there a short way to define BeamEnergy(Beam1_e, Beam2_e, particle1, particle2) class to incorporate different beam energies?

Thanks, Michael

[jncots]

Hello michael-pitt,

Do I understand correctly that you want, for example, collide proton with pz=1 TeV and Oxygen with pz=10 GeV? If so, you can do something like this:

import chromo
from chromo.kinematics import EventKinematics, EventFrame
from chromo.models import Sibyll23d
from chromo.constants import TeV, GeV

ekin = EventKinematics("p", "N", beam=(1*TeV, -10*GeV), frame=EventFrame.FIXED_TARGET)
gen = Sibyll23d(ekin)

for event in gen(100):
    event = event.final_state()
    print(event.en)
    ...

Note, that energy of final particles will be given in "EventFrame.FIXED_TARGET", where initial nucleus is at rest. You can choose EventFrame.CENTER_OF_MASS - center of mass of projectile and nucleon of nucleus. We haven't yet implemented transformation back to the initial frame "EventFrame.GENERIC". We can do this if needed.

[michael-pitt]

Hi @jncots

Thank you for your answer. The need is to run the simulations for collider experiments (like LHC) where both beam particles are accelerated. In case of pO collisions, I wish to have the following setup:

ekin = EventKinematics("p", "O", beam=(6.8*TeV, -8*6.8*GeV), frame=EventFrame.LAB_FRAME)

I'm not entirely sure, but is the current solution should be:

import chromo
import numpy as np
from matplotlib import pyplot as plt
from chromo.kinematics import EventKinematics, EventFrame
from chromo.constants import GeV, TeV
from particle import literals as lp
from pylorentz import Momentum4
import boost_histogram as bh

# the output frame is the proton rest frame so one need to boost it
boost = Momentum4(np.sqrt(lp.proton.mass**2 + (6.8*1e6)**2), 0, 0, 6.8*1e6) # note the different units (MeV <-> GeV)
ekin = EventKinematics("p", "O", beam=(6.8*TeV, -8*6.8*GeV), frame=EventFrame.FIXED_TARGET)
gen = chromo.models.Sibyll23d(ekin)

# book histograms
heta_fixed = bh.Histogram(bh.axis.Regular(50, -20, 20))
heta_boosted = bh.Histogram(bh.axis.Regular(50, -20, 20))

# simulate collisions and get rapidity
for event in gen(10):
    for par in event.final_state():
        par = Momentum4.m_eta_phi_pt(par.m, par.eta, par.phi, par.pt)
        heta_fixed.fill(par.eta)
        par = par.boost_particle(-boost) # proton at test, ion in positive direction
        heta_boosted.fill(par.eta)

# plot eta distribution of stable particles
fig, ax = plt.subplots(1, 2, figsize=(10, 5))
ax[0].bar(heta_fixed.axes[0].centers, heta_fixed.values(), width=heta_fixed.axes[0].widths);
ax[1].bar(heta_boosted.axes[0].centers, heta_boosted.values(), width=heta_boosted.axes[0].widths);

[afedynitch]

The correct way to simulate colliders is to provide the nucleon-nucleon center of mass energy, namely (6.8+Z/A*6.8)TeV, and get the output in center of mass frame.

[Xiaozizhuo]

你好

谢谢你的回答。需要为对撞机实验（如LHC）运行模拟，其中两个光束粒子都被加速。如果发生 pO 冲突，我希望进行以下设置：

ekin = EventKinematics("p", "O", beam=(6.8*TeV, -8*6.8*GeV), frame=EventFrame.LAB_FRAME)

我不完全确定，但当前的解决方案应该是：

import chromo
import numpy as np
from matplotlib import pyplot as plt
from chromo.kinematics import EventKinematics, EventFrame
from chromo.constants import GeV, TeV
from particle import literals as lp
from pylorentz import Momentum4
import boost_histogram as bh

# the output frame is the proton rest frame so one need to boost it
boost = Momentum4(np.sqrt(lp.proton.mass**2 + (6.8*1e6)**2), 0, 0, 6.8*1e6) # note the different units (MeV <-> GeV)
ekin = EventKinematics("p", "O", beam=(6.8*TeV, -8*6.8*GeV), frame=EventFrame.FIXED_TARGET)
gen = chromo.models.Sibyll23d(ekin)

# book histograms
heta_fixed = bh.Histogram(bh.axis.Regular(50, -20, 20))
heta_boosted = bh.Histogram(bh.axis.Regular(50, -20, 20))

# simulate collisions and get rapidity
for event in gen(10):
    for par in event.final_state():
        par = Momentum4.m_eta_phi_pt(par.m, par.eta, par.phi, par.pt)
        heta_fixed.fill(par.eta)
        par = par.boost_particle(-boost) # proton at test, ion in positive direction
        heta_boosted.fill(par.eta)

# plot eta distribution of stable particles
fig, ax = plt.subplots(1, 2, figsize=(10, 5))
ax[0].bar(heta_fixed.axes[0].centers, heta_fixed.values(), width=heta_fixed.axes[0].widths);
ax[1].bar(heta_boosted.axes[0].centers, heta_boosted.values(), width=heta_boosted.axes[0].widths);

Sorry to interrupt your discussion, I would like to ask a question. In your code, are these two modules in the package itself? pylorentz boost_histogram Or do you define it yourself? I don't find definitions for these two modules in my code file.

[jncots]

michael-pitt, I guess @afedynitch gave the right direction. If it is related to proton-ion collisions check https://physics.stackexchange.com/q/600966. In nucleon-nucleon CMS $E_{cms} \approx 2\sqrt{E_1E_2} \approx 2E_p\sqrt{\frac{Z_1Z_2}{A_1A2}}$, or for p-A: $E{cms}\approx 2E_p\sqrt{\frac{Z}{A}}$. In chromo:

e_beam = 6.8*TeV
emc = 2*e_beam*np.sqrt(8/16)
ekin = CenterOfMass(emc, "p", "O")

Using beams:

e_beam = 6.8*TeV
pz1 = e_beam
pz2 = -e_beam*8/16 # momentum per nucleon
ekin = EventKinematics("p", "O", beam=(pz1, pz2), frame=EventFrame.FIXED_TARGET)
gen = chromo.models.Sibyll23d(ekin)

If you want to have results in the rest frame of proton then following transformation should be used:

mnuc = (lp.proton.mass + lp.neutron.mass)/2*MeV
mprot = lp.proton.mass*MeV
inv_p1 = Momentum4(np.sqrt(mprot**2+pz1**2), 0, 0, -pz1) # inverse transform
inv_p2 = Momentum4(np.sqrt(mnuc**2+pz2**2), 0, 0, pz2)
transform_boost = inv_p2.boost_particle(inv_p1) # momentum inv_p2 in the rest frame of inv_p1

# book histograms
heta_fixed = bh.Histogram(bh.axis.Regular(50, -20, 20))
heta_boosted = bh.Histogram(bh.axis.Regular(50, -20, 20))

# simulate collisions and get rapidity

for event in gen(100):
    evt_fin = event.final_state()
    heta_fixed.fill(evt_fin.eta) # much faster than iteration
    for par in evt_fin: # not very efficient
        par = Momentum4(par.en, par.px, par.py, par.pz)
       # Transform from nucleus frame to initial frame, and from initial frame to proton frame
        par = par.boost_particle(transform_boost) # proton at rest, ion in negative direction
        heta_boosted.fill(par.eta)

# plot eta distribution of stable particles
fig, ax = plt.subplots(1, 2, figsize=(10, 5))
ax[0].bar(heta_fixed.axes[0].centers, heta_fixed.values(), width=heta_fixed.axes[0].widths);
ax[1].bar(heta_boosted.axes[0].centers, heta_boosted.values(), width=heta_boosted.axes[0].widths);

[jncots]

@Xiaozizhuo, pylorentz, boost_histogram can be installed as usual:

pip install pylorentz boost_histogram

They are not a part of Chromo.

By the way, official release of chromo 0.4 is on PyPI now. You can install/update it as:

pip uninstall chromo
pip install chromo

[Xiaozizhuo]

Thank you for your help, it will be of great use to my later studies

[michael-pitt]

Hi @afedynitch,

Thanks for the info. It is great to know that the beam parameters are nucleon energy (similar to many other MC generators). The problem with the center of mass frame is that it is not the lab frame. For NN collisions with beam=(6.8*TeV, -0.5*6.8*TeV) in the center of the mass frame, both nucleons will have 5.1*TeV, so I need to boost the particles to the lab frame as far as I understand.

It will be good to add an example of running pA collisions, as many are using these models, and the EventFrame.GENERIC option will solve the CM - lab frame issue.

@jncots @afedynitch I noticed something strange: chromo.kinematics.CenterOfMass( (6.8 + 8/16*6.8)*TeV, "p", "O") and chromo.kinematics.EventKinematics("p", "O", beam=(6.8*TeV, -0.5*6.8*TeV), frame=EventFrame.CENTER_OF_MASS ) gives different results (different beam energies), I don't know what I'm doing not correctly.

@Xiaozizhuo, I managed to avoid using pylorentz and boost_histogram, and switched using vector. Here is my updated code snippet

import chromo
import numpy as np
from matplotlib import pyplot as plt
from chromo.kinematics import EventKinematics, EventFrame
from chromo.constants import GeV, TeV
from tqdm import tqdm
import vector

# The output frame is the CM rest frame, so one needs to boost it to the lab frame
betaZ = ( 6.8**2 - 5.1**2) / ( 6.8**2 + 5.1**2)  # 5.1 = (6.8 + 8/16*6.8) / 2
boost = vector.obj(x=0, y=0, z=betaZ)

ekin = chromo.kinematics.CenterOfMass( (6.8 + 8/16*6.8)*TeV, "p", "O")
#ekin = EventKinematics("p", "O", beam=(6.8*TeV, -0.5*6.8*TeV) # not working

gen = chromo.models.Sibyll23d(ekin)
n_events = 1000

# create arrays
eta_fixed = [np.array([], dtype=float)]*3
eta_boosted = [np.array([], dtype=float)]*3

# simulate collisions and get rapidity
taga=False
for event in tqdm(gen(n_events), total=n_events):
    for par in event.final_state():
        pid = np.abs(par.pid)
        if pid in [130, 310, 311, 321]: pid=0
        elif pid in [211, 111]: pid=1
        elif pid in [2212]: pid=2
        else: continue
        par = vector.obj(pt=par.pt, phi=par.phi, eta=par.eta, mass=par.m)
        eta_fixed[pid] = np.append(eta_fixed[pid],par.eta)
        par = par.boost(boost) # boost to the lab frame
        eta_boosted[pid] = np.append(eta_boosted[pid], par.eta)

# plot eta distribution of stable particles
fig, ax = plt.subplots(1, 2, figsize=(10, 5))
for i, pname in enumerate(['K',r'$\pi$','p']):
  ax[0].hist(eta_fixed[i], bins=np.linspace(-15,15,50), label=pname, alpha = 0.5, density=True); ax[0].set_xlabel('particle eta'); ax[0].set_title(r'CM frame: $p\rightarrow$ $\leftarrow O$');
  ax[1].hist(eta_boosted[i], bins=np.linspace(-15,15,50), label=pname, alpha = 0.5, density=True); ax[1].set_xlabel('particle eta'); ax[1].set_title(r'Lab frame: $p\longrightarrow$ $\leftarrow O$');
ax[0].legend()
plt.savefig('set_pA.png')

[jncots]

I noticed something strange ... gives different results

These are approximately equivalent settings:

ekin1 = chromo.kinematics.CenterOfMass( 2*6.8*np.sqrt(0.5)*TeV, "p", "O")
# and 
ekin2 = chromo.kinematics.EventKinematics("p", "O", beam=(6.8*TeV, -0.5*6.8*TeV), frame=EventFrame.CENTER_OF_MASS )

print(f"rel_error_elab = {(ekin1.elab-ekin2.elab)/ekin2.elab},\nrel_error_ecm = {(ekin1.ecm-ekin2.ecm)/ekin2.ecm}")

Note that $E_{cms}\approx 2Ep\sqrt{\frac{Z}{A}}$. This is approximate expression, for $E \gg m$. See previous comment. Exact expression is also easy to derive: $E{cms} = \sqrt{m^2_1 + m^2_2 + 2E_1E_2\left(1+\sqrt{1-\frac{1}{\gamma^2_1}}\sqrt{1-\frac{1}{\gamma^2_2}}\right)}\approx\sqrt{m^2_1 + m^2_2 + 4E_1E_2}$

[jncots]

It will be good to add an example of running pA collisions, as many are using these models, and the EventFrame.GENERIC option will solve the CM - lab frame issue.

Ok, we will add generic transformations.

[afedynitch]

@michael-pitt, it depends what you aim to do. If you are setting up a low level simulation of an experiment, then maybe what you propose is the correct thing to do. If you just want to compare to published results, then the published data is already boosted back to the CMS. Note that fixed target data, from e.g. NA61 is also published in CMS.

Nonetheless, as @jncots mentioned above, we should implement generic transformations such that the initial conditions can be defined by 2 4vectors.

[michael-pitt]

Hi @jncots, @afedynitch

Thank a lot for the feedback, and #184

@michael-pitt, I guess @afedynitch gave the right direction. If it is related to proton-ion collisions check https://physics.stackexchange.com/q/600966. In nucleon-nucleon CMS Ecms≈2E1E2≈2EpZ1Z2A1A2, or for p-A: Ecms≈2EpZA. In chromo:

e_beam = 6.8*TeV
emc = 2*e_beam*np.sqrt(8/16)
ekin = CenterOfMass(emc, "p", "O")

Oh, I missed your previous comment. You are right, ECMS is not a sum of nucleon energies in case of ions but $E_{CMS}=2\sqrt{E_aE_b}$ where $E_a$ and $E_b$ is the energy per nucleon. This could be another argument for providing an example in chromo/examples, :satisfied:

@michael-pitt , it depends what you aim to do. If you are setting up a low level simulation of an experiment, then maybe what you propose is the correct thing to do. If you just want to compare to published results, then the published data is already boosted back to the CMS. Note that fixed target data, from e.g. NA61 is also published in CMS.

I'm trying to produce events from a collider experiment in a detector rest frame. This is common practice since particles are processed later through a detector simulation, so all particles are fixed to the lab frame.

I've copied bellow a corrected version of the code (to avoid confusion with the previous incorrect implementation)

import chromo
import numpy as np
from matplotlib import pyplot as plt
from chromo.kinematics import EventKinematics, EventFrame
from chromo.constants import GeV, TeV
from tqdm import tqdm
import boost_histogram as bh
import vector

# define collision energy in center of mass frame
e_beam = 6.8*TeV
pz1 = e_beam
pz2 = -e_beam*8./16. # momentum per nucleon
ekin = EventKinematics("p", "O", beam=(pz1, pz2), frame=EventFrame.CENTER_OF_MASS )
#ekin = chromo.kinematics.CenterOfMass(2*e_beam*np.sqrt(8/16), "p", "O")

# define the boost from CM frame to lab frame:
betaZ = ( 1 - 8./16. ) / ( 1 + 8./16. )
boost = vector.obj(x=0, y=0, z=betaZ)

gen = chromo.models.Sibyll23d(ekin)
n_events = 10000

# book histograms
heta_cms = [bh.Histogram(bh.axis.Regular(50, -15, 15)), bh.Histogram(bh.axis.Regular(50, -15, 15)), bh.Histogram(bh.axis.Regular(50, -15, 15))]
heta_lab = [bh.Histogram(bh.axis.Regular(50, -15, 15)), bh.Histogram(bh.axis.Regular(50, -15, 15)), bh.Histogram(bh.axis.Regular(50, -15, 15))]

# simulate collisions and get rapidity
for event in tqdm(gen(n_events), total=n_events):
    for par in event.final_state():
        pid = np.abs(par.pid)
        if pid in [130, 310, 311, 321]: pid=0
        elif pid in [211, 111]: pid=1
        elif pid in [2212]: pid=2;
        else: continue
        par = vector.obj(pt=par.pt, phi=par.phi, eta=par.eta, mass=par.m)
        heta_cms[pid].fill(par.rapidity)
        par = par.boost(boost) # boost to lab frame
        heta_lab[pid].fill(par.rapidity)

# plot eta distribution of stable particles
fig, ax = plt.subplots(1, 2, figsize=(10, 5), sharey=True)
for i, pname in enumerate(['K',r'$\pi$','p']):
  ax[0].bar(heta_cms[i].axes[0].centers, heta_cms[i].values()/heta_cms[i].sum(), width=heta_cms[i].axes[0].widths, label=pname, alpha = 0.5); 
  ax[0].set_ylabel('normalized'); ax[0].set_xlabel('particle rapidity'); ax[0].set_title(r'CM frame: 4.8 TeV p$\rightarrow$ $\leftarrow$ 76.8 TeV O');
  ax[1].bar(heta_lab[i].axes[0].centers, heta_lab[i].values()/heta_lab[i].sum(), width=heta_lab[i].axes[0].widths, alpha = 0.5); ax[1].set_xlabel('particle rapidity'); ax[1].set_title(r'Lab frame: 6.8 TeV p$\rightarrow$ $\leftarrow$ 64.4 TeV O');
ax[0].legend()
plt.savefig('set_pA.png')