Solving Lorentz equation and the effect of a radial field in HELIOS

Leave a comment

Haven’t posted anything for a while, was working on extracting carbon isotopes ESPS, the FSU DAQ, and SOLARIS DAQ for the NUMEN collaboration, went to Hawaii for DNP 2023, and got sick for a week…. Anyway, a friend was having trouble doing a HELIOS simulation with a non-uniform field in GEANT4, and I am curious what is the effect of a non-uniform field, so, I decided to do my own simulation. The core of the simulation is solving the Lorentz equation \vec{F} = Qe \vec{v} \times \vec{B} . This is a 2nd order coupled equation. Well, I already solved it before, so it would be not a huge work.


In this post, I already show that using RK4 method to solve the Lorentz equation. Here, I show the method using the Scipy library in Python and also NDSolve[] in Mathematica.

#Define magnetic field
def Bx(x, y, z):
  return  x/ (1 + math.exp( x**2 + y**2) )
def By(x, y, z):
  return  y/ (1 + math.exp( x**2 + y**2) )
def Bz(x, y, z):
  return -2.5 / (1 + math.exp((math.sqrt(x**2 + y**2) - 0.5)/0.05)) / (1 + math.exp((z**2 - 1**2)/0.1))

from scipy.integrate import odeint

def FieldEquations(w, t, p):
  x , y, z, vx, vy, vz = w

  f = [vx,
       vy,
       vz,
       p * ( vy * Bz(x,y,z) - vz * By(x,y,z)),
       p * ( vz * Bx(x,y,z) - vx * Bz(x,y,z)),
       p * ( vx * By(x,y,z) - vy * Bx(x,y,z))]
  return f

def SolveTrace( helios, thetaCM_deg, phiCM_deg, r0_meter, tMax_ns):
  vx0, vy0, vz0 = helios.velocity_b(thetaCM_deg, phiCM_deg)
  w0 = [r0_meter[0], r0_meter[1], r0_meter[2], vx0, vy0, vz0]

  # ODE solver parameters
  abserr = 1.0e-8
  relerr = 1.0e-6
  stoptime = tMax_ns
  numpoints = 500

  # Create the time samples for the output of the ODE solver.
  t = [stoptime * float(i) / (numpoints - 1) for i in range(numpoints)]

  wsol = odeint(FieldEquations, w0, t, args=(helios.QMb,), atol=abserr, rtol=relerr)

  return [t, wsol]

The code is very simple, we first define the magnetic field components, next, in the FieldEquation, the w is the coordinate vector, and f is the corresponding equation for the coordinate vector, for the x-component, which is

\displaystyle \frac{dx}{dt} = v_x \\ \frac{dv_x}{dt} = p ( v_y B_z - v_z B_y )

Similar for the y- and z-components. The variable p is the charge-to-mass ratio. The SolveTrace will use the helios class, which is not shown here, but it basically gives the velocity vector for the initial velocity. The output is the list of the time and the wsol. the wsol is an array of array of x, y, z, vx, vy, and vz.

In Mathematica,

fieldX[x_, y_, z_] := 0.1 If[x == 0, 1, x/Sqrt[x^2 + y^2] Exp[-2 (x^2 + y^2)]]
fieldY[x_, y_, z_] := 0.1 If[y == 0, 1, y/Sqrt[x^2 + y^2] Exp[-2 (x^2 + y^2)]]
fieldZ[x_, y_, z_] := -2.5

G2x[t_, r_, v_] := v[[2]] fieldZ[r[[1]], r[[2]], r[[3]]] - v[[3]] fieldY[r[[1]], r[[2]], r[[3]]];
G2y[t_, r_, v_] := v[[3]] fieldX[r[[1]], r[[2]], r[[3]]] - v[[1]] fieldZ[r[[1]], r[[2]], r[[3]]];
G2z[t_, r_, v_] := v[[1]] fieldY[r[[1]], r[[2]], r[[3]]] - v[[2]] fieldX[r[[1]], r[[2]], r[[3]]]

SolvePath[ChargetoMass_, G2x_Symbol, G2y_Symbol, G2z_Symbol,  r0_, v0_, tMax_] :=
 NDSolve[{
   x''[t] == ChargetoMass G2x[t, { x[t], y[t], z[t]}, { x'[t], y'[t], z'[t]}],
   y''[t] == ChargetoMass G2y[t, { x[t], y[t], z[t]}, { x'[t], y'[t], z'[t]}],
   z''[t] == ChargetoMass G2z[t, { x[t], y[t], z[t]}, { x'[t], y'[t], z'[t]}],
   x[0] == r0[[1]], x'[0] == v0[[1]],
   y[0] == r0[[2]], y'[0] == v0[[2]],
   z[0] == r0[[3]], z'[0] == v0[[3]]
   }, {x[t], y[t], z[t]}, {t, 0, tMax}]

The code defines the field components first, the G2x is the (\vec{v} \times \vec{B} )_x . After that, use the NDSolve to solve for x(t), y(t), z(t) .

A small note on the unit. We are going to solve it in units of meters and nano-seconds. So, the initial position and velocity are in meter and meter/nano-sec respectively. and the charge-to-mass ratio is in SI unit, so, using Coulomb over kg. The charge-to-mass ratio also needs to be multiplied by 10^{-9} because the left-hand side is 2nd derivative of time in nano-sec, while the right-hand side needs to be multiplied by 10^{-9} to give the correct scale. (think about the original Lorentz equation is in mater and sec)


Now, we are in a position to solve the locus of the charged particle for any field. In particular, I like to find out the effect of a radial component.

I set the radial field to be B_\rho (\rho, z) = k^2 \exp(-8 \rho^2 ) , where k is the strength of the radial component. The radial field is strong at small \rho and decay with \rho.

\displaystyle B_x(x, y, z) = k \frac{x}{\sqrt{x^2+y^2}} \exp(-4 (x^2+y^2)) \\ B_y(x, y, z) = k \frac{y}{\sqrt{x^2+y^2}} \exp(-4 (x^2+y^2)) \\ B_z (x,y, z) = - 2.5

For the reaction, I pick 15N(d,p) at 10 MeV/u. At \theta_{CM} = 20 ^\circ, The velocity vector is ( 0.0195276, 0, -0.0147025 ) [meter/ns]. The locus when k = 0.1 is

Interestingly, the particle experiences a positive z-acceleration. Imagine the particle only moves in the x-y plane, with a radial component of the field, the particle experiences a positive z-force.

Another effect is the distortion of the x-y projection. Without the radial field, the x-y projection should be a perfect circle. Since the radial field gives a +z-force, and the z-field gives a radially inward force. These two forces combined will twist the path.

Here is the locus for k = 0.5

as expected, a stronger radial field, the stronger positive z-force, and the particle change the moving direction more quickly. Also, the x-y project is more distorted.

Compare Non-relativistic and relativistic calculation in HELIOS

Leave a comment

The classical result for inverse kinematic is, where m is in SI unit,

\displaystyle z = T_c ( -u \cos\theta + V_{cm} ) \\ E = \frac{1}{2} m(u^2 + V_{cm}^2) - m u V_{cm} \cos\theta

The relativistic result is where energy is in MeV unit,

\displaystyle z = \frac{2\pi}{cqB} \left( \gamma \beta E_{cm} - \gamma p \cos\theta\right) \\ E = \gamma E_{cm} - \gamma \beta p \cos\theta

The classical limit can be reached as

\displaystyle \gamma \rightarrow 1, ~~ \beta \rightarrow \frac{V_{cm}}{c}

and the classical cyclotron period is T_c = 2\pi m / qB . We have:

\displaystyle z = \frac{2\pi}{qB} \left( V_{cm} \frac{E_{cm}}{c^2} -  \frac{p}{c} \cos\theta\right) \\ E =  E_{cm} - V_{cm} \frac{p}{c} \cos\theta

We make more approximations,

\displaystyle E_{cm} \approx mc^2, ~~ p \approx mu c

Then

\displaystyle z = \frac{2\pi m}{qB} \left( V_{cm} - u\cos\theta\right) \\ E =  E_{cm} - V_{cm} mu \cos\theta

We kept the E_{cm} in the 2nd equation because we have to subtract mc^2 at the end.

\displaystyle z = T_c \left( V_{cm} - u\cos\theta\right) \\ E =  E_{cm} - V_{cm} mu \cos\theta

We can see that, the z-position assumes the form of the classical case, and by comparison,

\displaystyle E_{cm} \approx mc^2 + \frac{1}{2} m(u^2 + V_{cm}^2)

That means, in the classical limit, compared to relativistic, the transformation is simply adding CM frame motion.

The intermediate step in eliminating the \cos\theta from either case is

\displaystyle E = \frac{1}{2} m(u^2 + V_{cm}^2) - mV_{cm}^2 + \frac{m V_{cm}}{Tc} z\\

\displaystyle E = \frac{1}{2} m(u^2 - V_{cm}^2) + \frac{m V_{cm}}{T_c} z\\

Compared to the final result from a fully relativistic

\displaystyle E = \frac{1}{\gamma} E_{cm} + \frac{mc\beta}{T_c} z , ~~ T_c = \frac{2\pi}{B} \frac{m}{q}

It may be misleading that E_{cm} \approx mc^2 + \frac{1}{2} m(u^2 - V_{cm}^2) .

HELIOS Kinematic in Non-relativistic limit

Leave a comment

The reaction is 2(1, a) b. where particle 1 with kinetic energy of T, velocity of v_1 = \sqrt{2T/m_1} hit particle 2, resulting particle a and b, such that 1 becomes a, 2 becomes b. The velocity of the CM frame is

\displaystyle V_{cm} = \frac{m_1}{m_1+m_2} v_1

In the CM frame, the particle velocities are

\displaystyle u_1 = v_1 - V_{cm} = \frac{m_2}{m_1+m_2} v_1 = \frac{m_2}{m_1}V_{cm}

\displaystyle u_2 = 0 - V_{cm} = - \frac{m_1}{m_1+m_2} v_1 = - V_{cm}

The total KE in CM frame is

\displaystyle T_{cm} = \frac{1}{2}m_1 u_1^2 + \frac{1}{2}m_2 u_2^2 = \frac{1}{2}\frac{m_1m_2}{m_1+m_2} v_1^2 = \frac{1}{2} \frac{m_2}{m_1} (m_1+m_2) V_{cm}^2

The T_{cm} can be expressed as T as

\displaystyle T_{cm} = \frac{m_2}{m_1+m_2} T

After the scattering, the kinematics of the particle a and b must conserve the energy and momentum, i.e.

\displaystyle  \frac{1}{2}m_a u_a^2 + \frac{1}{2}m_b u_b^2 = T_{cm} + Q

\displaystyle m_a u_a + m_b u_b = 0

\displaystyle Q = m_1 + m_2 - m_a - m_b

The solution for u_a is

\displaystyle u_a = \sqrt{\frac{2m_b Q}{m_a(m_a+m_b)} + M V_{cm}^2}, ~~ M = \frac{m_2}{m_1} \frac{m_1+m_2}{m_a+m_b}\frac{m_b}{m_a}

Move to Lab frame,

\displaystyle v_a = \begin{pmatrix} u_a \cos\theta + V_{cm} \\ u_a \sin\theta \end{pmatrix}

The particle a carry charge Q and moving in a magnetic field of strength B in the z-direction. The distance in the beam direction after time T_c is

\displaystyle \frac{z}{T_c}  = u_a \cos \theta  + V_{cm}

The KE of the particle a in Lab frame is

\displaystyle E_a = \frac{1}{2} m_a \left( u_a^2  + V_{cm}^2 + 2 u_a V_{cm} \cos\theta\right)

eliminate the center of mass scattering angle using the z-position, we have

\displaystyle E_a = T_a  - \frac{1}{2} m_a V_{cm}^2 +  \frac{m_a V_{cm}}{T_c} z, ~~~ T_a = \frac{1}{2} m_a u_a^2

The T_a is the KE of the particle a in CM frame. The excitation of particle b is implicitly included in the term T_a through the Q value. The 2nd term depends on the initial condition, and the KE of particle a is linear of the z position.

Once the term u_a deduced, the Q value is

\displaystyle Q = \frac{m_a+m_b}{m_b}T_a - \frac{1}{2} \frac{m_2}{m_1} (m_1+m_2) V_{cm}^2 =  \frac{m_a+m_b}{m_b}T_a - T_{cm}


In above calculation, I assumed the particle 1 (or a) is lighter then particle 2 (or b). so, after scattering, the particle a is moving “forward”.

In inverse kinematics, the light particle is b. so

\displaystyle u_b =\frac{m_a}{m_b} u_a = \frac{m_a}{m_b} \sqrt{ \frac{2 m_b Q}{ m_a(m_a+m_b)} +  M V_{cm}^2}

in the lab frame,

\displaystyle v_b = \begin{pmatrix} u_b \cos\theta + V_{cm} \\ u_b \sin\theta \end{pmatrix}

The KE of particle b is

\displaystyle E_b = \frac{1}{2} m_b \left( u_b^2  + V_{cm}^2 - 2 m_b V_{cm} \cos\theta\right)

And the E-z equation is

\displaystyle E_b = T_b - \frac{1}{2} m_b V_{cm}^2 + \frac{m_b V_{cm}}{T_c} z, ~~~~T_b =  \frac{1}{2} m_b u_b^2

The Q value is

\displaystyle Q = \frac{m_a+m_b}{m_b} T_b - T_{cm}


In fact, from the balance of energy and the momentum

\displaystyle T_a + T_b = Q + T_{cm} , ~~ m_a u_a = m_b u_b

write the momentum squared as

\displaystyle m_a^2 u_a^2 = 2 m_a T_a = 2 m_b T_b

By eliminating T_a or T_b, We have either

\displaystyle T_a + \frac{m_a}{m_b}T_a = Q + T_{cm}

or

\displaystyle \frac{m_b}{m_a}T_b + T_b = Q + T_{cm}

Q-value and z-position in HELIOS

Leave a comment

in HELIOS, I always have a feeling that the Q-value is related to the z-position. Look at the Q-value,

\displaystyle Q = m_A + m_a - m_b - m_B - E_x,

where E_x is the excitation energy of the heavy recoil B. And we know that the higher the excitation energy, the z-position will shift toward the positive z-position. But that is for a given reaction. Is there any general or universal relationship between the Q-value and the z-position?

Instead of solving the equations and finding out the relation, I use simulation, put various reactions, and beam nuclei (6,9Li, 10B, 12C, 16N, 16,20O, 19F, 25Mg, 30,38Si, 40,48Ca, 96Zr, 158Tb, 208Pb, 238U ) adjusting the E_x from 0 to 2 MeV to mimic the change of Q-value. And I fixed the z-position at \theta_{cm} = 10 \circ , the beam energy of 10 MeV/u, and the magnetic field was 2.5 T. Here is the result:

For the given reaction, say the two neutron transfer reactions (d,p) and (t,d) has a similar trend. And (d,3He) and (d,a) are also on the different sections of the same curve. In general, all curve, except the (d,t), has a similar slope, roughly -32 mm/MeV, the (d,p) and (t,d) is a bit steeper with a slope of -37 mm/MeV.

Also, in the 2-neutron-adding (t,p) reaction, the z-position (-800 mm) at Q =0 is roughly double of the 1-neutron-adding (d,p) or (t,d) reaction, which is ( ~500 mm) . And charge exchange reaction (t, 3He) is almost at z = 0. The proton-removal reactions (d, 3He) and (t, a) share a similar verticle offset of 400 mm. and neutron-removal and pn-removal reactions have z-position at ~ 700 mm. In general, adding reactions has negative z, and removal reaction has positive z for zero Q-value.

Last note, with the higher magnetic field, the curves will move toward z = 0, and the slope will also change, for example, the (d,p) reaction

We can see the effect of the B-field is stronger for a higher Q-value reaction for (d,p) reaction. In general, we can foresee there is a Q-value that is unaffected by the B-field, and any Q-value further away from this special Q-value will be affected by the B-field the more. In fact, since the “contraction” is around z=0, the special Q-value that is unaffected by the B-field is near the curve at z=0. For (d,p) and (t,d) reactions, this special Q-value is around -9 MeV. For (t, 3He) reaction, this special Q-value is ~ 0 MeV. For (t, a) reaction, it is around 11 MeV.

And If I change the beam energy but keep the B-field to 2.5 T for the (d,p) reaction:

Higher the beam higher, the higher the energy of the light recoil, and the z-positon move further away.

Web-based HELIOSmatics

Leave a comment

There are 2 main tools for the kinematic simulation or calculation for the HELIOS spectrometer. One is using MS Excel, and another one is using CERN ROOT with Monte Carlo events generation.

They work great after downloaded and with MS Excel, not Libre spreadsheet. But a problem is, whenever the code got updated, the user need to update it or download it again. Another disadvantage is on the CERN ROOT simulation, it requires few configuration files. When users need to preform multiple difference reactions and settings, it is not so convivence. In fact, the fundamental layer of the simulation code support different file names, but the GUI does not.

Now I moved the simulation on the web. The calculation mostly done with JavaScript, and plotting by Plotly. The DWBA and Monte Carlo simulation is using CGI python interface to reuse the code developed in C++, like retrieving the mass, calling the DWBA code, etc.

The reason it can be put on the web is the collaboration between FSU and FRIB on the SOLARIS project.

Mathematics of Binning

Leave a comment

Given that a distribution f(x), it populate events on the demain from x_1 to x_2. Let say, we like to plot the histogram with binning of \Delta. Mathematically, we are doing an integration

\displaystyle b_i = b(x_1+i \Delta)  = \int_{x_1 + i \Delta}^{x_1 + (i+1)\Delta} f(x) dx = F\left(x_1 + (i+1) \Delta \right) - F(x_1 + i \Delta)

If we divide each bin with \Delta

\displaystyle \frac{b_i}{\Delta} =  \frac{1}{\Delta}  \left( F\left(x_1 + (i+1) \Delta \right) - F(x_1 + i \Delta)  \right) \approx f(x_1 + i \Delta)

or

\displaystyle b(x) = \Delta f(x)

That means, when we “translate” the histogram into a scattering plot and we can treat it as a sampling of the distribution. In other word, depends on how the y-axis being labeled, the same data represents different meanings.


In nuclear physics, the differential cross section is the intrinsic quantity that “generates” the histogram. For example, a measured angular distribution is

\displaystyle  f(\theta) = \frac{d\sigma}{d\Omega},   b_i \propto \int_{\Omega_i} f(\theta) \sin(\theta) d\theta d\phi

Usually, the label of the y-axis of the histogram is “count”, “count in angle bin”. Or, less accurate to be “count per angle”, or “count / \Omega “. If we treat the histogram as a scattering plot,

\displaystyle b(\theta) \approx f(\theta) \sin(\theta)

and then normalized the height with \sin(\theta), the beam intensity, and target thickness etc, then we can label the y-axis as \frac{d\sigma}{d\Omega} \textrm{[mb/sr]} .


In Helios, things get a bit tricky. Because of the magnetic field, we have

\displaystyle z \propto \cos(\theta)  \rightarrow dz = \sin(\theta) d\theta

Thus, if we treat the count for a detector, which integrate the differential cross section over the acceptance of a single detector,

\displaystyle b_i \propto \int_{\Omega_i} \frac{d\sigma}{d\Omega} d\cos(\theta) d\phi = \int_{\Omega_i} \frac{d\sigma}{d\Omega} dz d\phi

Following the above treatment, if we treat the histogram as a scattering plot, thus, we DONNOT need to normalize the data with \sin(\theta) , and get the differential cross section.

Finite detector radius II

Leave a comment

What had I done in last 2 months!?? Modified/improved the RAISOR DAQ for MCP, prepared and running 29Al(d,p). 16N isomer measurement (early Oct), target degradation data, DNP2020 talk, ISOLDE talk, review the HELIOS calculation and using Mathematica to visualize it, help MUSIC to setup the digitizer DAQ (kind of fail… ), that’s all? need to work harder…..

Calculation of the 29Al(d,p) @ 10.45 MeV.u, 2.5 T. the recoil XY scaled up 10 times.

In this post, we addressed the finite detector radius correction. When we plotting the e-z plot, we assumed \Delta\phi = \phi-\phi_d = \pi as the orbit radius \rho >> a, where a is the detector radius, \phi is the scattering azimuthal angle of the light particle, and \phi_D is the detector plane angle.

As we can expect, when \theta_{cm} < 10 ^{\circ} and small, the orbit radius \rho is also small, and the above approximation is far from correct. When \theta_{cm} is small, we can imagine the \Delta \phi = \phi - \phi_d < \pi , for example, as shown in below,

Lets recall the finite-detector size correction

\displaystyle z_{Hit} = \frac{z_0(\theta_{cm})}{2\pi} \left( \Delta\phi +\pi - \sin^{-1} \left( \frac{a}{\rho(\theta_{cm})} - \sin(\Delta \phi)  \right)  \right)

If we treat \Delta\phi and \theta_{cm} as variable, we have following e-z plot,

We can see that, there is a \Delta \phi_X value for a “extreme” z-position, thus, it is natural to differentiate the z_{Hit} with respect to \Delta \phi, and solve for \Delta \phi_X . The result is,

\displaystyle \sin (\Delta \phi_X) = \pm \frac{a}{2\rho}

In here, the \pm sign is for magnetic field parallels or anti-parallels the beam direction. Since we know that \Delta \phi_X must be close to \pi , we have the solution

\displaystyle \Delta \phi_X = \pi \mp \sin^{-1}\left( \frac{a}{2\rho(\theta_{cm})} \right)

From the solution, we can see that when a \sim \rho , the \Delta \phi_X different from pi by a lot. Substitute back the the z-position, for n = 1 or only 1-loop

\displaystyle z_X = z_0 \left( 1 - \frac{1}{\pi} \sin^{-1}\left( \frac{a}{2\rho(\theta_{cm})} \right) \right)

Compare to the \Delta \phi = \pi approximation,

\displaystyle z_\pi = z_0 \left( 1 - \frac{1}{2\pi} \sin^{-1}\left( \frac{a}{\rho(\theta_{cm})} \right) \right)

The only difference is the 2 inside or outside the arcsin. In the following plot, the configuration is the same, and we added the red curve for the z_X .

The effect is more significant when a is large, for example, a = 50 mm. In the following plot, the blue curve is z_\pi , the orange curve is z_{140} , the red is z_X, and the green curve is out of range.

We also did a simulation for a = 50 mm, 4-side. In the following plot shows the result. The z_X (red) and z_\pi (green) are plotted as well. We can see the z_X line match the simulated data better.

It is worth to note that, the end of the z_\pi curve is at the edge of the e-z distribution. In fact, when we plot a contour of many \Delta\phi , it span the e-z distribution. For a 6-side array, the span will be smaller. If the detector is cylindrical, the span will be gone and reduces to the z_X. In that case, z_\pi is incorrect.

We also simulated that the excitation energy of the the heavy recoil follows normal distribution with mean of 0 MeV and sigma of 0.1 MeV. (there is a small deviation near the bending. )

Posters for Postdoc symposium

Leave a comment

I made two posters, one for the experiment on 207Hg and the other is about real-time beam diagnostic.

2019_postdoc_207Hg-1.png

2019_postdoc_HELIOS-1.png

Argonne HELIOS spectrometer & its scientific success

Leave a comment

I have a chance to give a small seminar in Hong Kong University. Here I attached the ppt as pdf file.

argonne helios spectrometer &amp; its scientific success

HKU_talk_cover.PNG

The talk is an overview of the HELIOS spectrometer in Argonne National Laboratory ( which I am mainly working on). This talk also help me what is done and what is not explored. Hope you will find it interesting.

Updated android program – Rel.Cal.

Leave a comment

A page for HELIOS kinematics was added in the previous program.

Click HERE to get the apl file

Here is the screenshot

screenshot_20190109-215434_rel cal

The plot has a lot to improve, say

  • auto-adjustment of the origin.
  • Added a dot to indicate the current thetaCM
  • able to change the line for constant thetaCM

Older Entries