Homework Solutions

parName <- c("electron", "muon", "proton", "deuteron")
mc2 <- c(0.510999,    105.658,  938.272,  1875.613)    # MeV
amm <- c(0.00115965, 0.00116592, 1.7928,  -0.14299)

Monday (due Tuesday)

  1. Magical Momenta. The anomalous magnetic moments of the electron, muon, proton, deuteron are \(a\) = (0.0011596, 0.0011659, 1.7928, -0.14299), respectively. Their rest mass energies are (0.510999, 105.658, 938.272, 1875.613) MeV, respectively. For each particle,

    1. determine its magic momentum and the corresponding kinetic energy.

    2. estimate the size of a storage ring operating at its magic momentum assuming a continuous magnetic field with strength 1.5 T.

bg = 1/sqrt(amm)
gam = sqrt(1+bg^2)
pmag = bg*mc2                 # MeV/c
KEmag = (gam-1)*mc2           # MeV
brho = pmag/1000*10/2.9979    # T-m
R0 = brho/1.5                 # m
C0 = 2*pi*R0                  # m
out = cbind.data.frame(parName,round(mc2,2),amm,round(pmag,2),round(KEmag,2),round(R0,3))
kable(out,col.names = c("Particle","mc2 (MeV)","a","pmag (MeV/c)","KEmag (MeV)","R0 (m)"))
Particle mc2 (MeV) a pmag (MeV/c) KEmag (MeV) R0 (m)
electron 0.51 0.0011596 15.01 14.50 0.033
muon 105.66 0.0011659 3094.34 2990.49 6.881
proton 938.27 1.7928000 700.75 232.80 1.558
deuteron 1875.61 -0.1429900 NaN NaN NaN
  1. Thomas Precession. The Thomas precession equation \[ \frac{d\bf s}{dt}=\frac{e}{mc}{\bf s\times}\left[\left(\frac{g}{2}-1+\frac{1}{\gamma}\right){\bf B} - \left(\frac{g}{2}-1\right)\frac{\gamma}{\gamma+1}({\bf\beta\cdot B)\beta}-\left(\frac{g}{2}-\frac{\gamma}{\gamma+1}\right){\rm \beta\times E} \right] \] relates the rate of change of the polarizaton \(({\bf s})\) to the electro-magnetic fields. When \({\bf B}\) is in the y-direction, \(\vec\beta = \beta {\bf \hat z}\), and \({\bf E}=0\), the polarization precesses about the y-axis with frequency \(\omega_p = \frac{e}{mc}\left(\frac{g}{2}-1+\frac{1}{\gamma}\right)B\). The cyclotron frequency \(\omega_c = \frac{eB}{mc\gamma}\) and \(\omega_a = \omega_p-\omega_c\).

    1. Suppose that the velocity is constant with direction \(\vec\beta = \beta (\sin\psi {\bf\hat y}+\cos\psi{\bf\hat z})\). Determine \(\vec\omega_p\), (magnitude and direction).

    2. What is the cyclotron frequency \(\vec{\bf \omega}_c\)?

    3. What are \({\bf\vec\omega_p}\) and \({\bf\vec\omega_c}\) in the limits where \(\psi=0\) and \(\psi=\pi/2\)?

  1. Betatron Tunes. Consider the \(g\)-2 storage ring, with ideal circular trajectory of radius \(R_0\) = 7.112 m.

    1. Given the magic momentum of the muon, compute the magnetic field required to acheive this radius.

    2. The electrostatic quadrupoles in the ring have electric fields that are zero on the design trajectory and, to lowest order, vary linearly with \(x\) and \(y\). The high voltage plates are located at \(\pm 50\) mm from the central trajectory. For a high voltage of \(\pm 20\) kV, compute the field gradient, \(E'\equiv \partial E_x/\partial x\) and find the “\(K\)” values of Hill’s Equation (\(K_x\) and \(K_y\)) for

      1. a region of length \(\ell_b\) with only \(B\) (no \(E\) field) and

      2. a region of length \(\ell_q\) with both \(B\) and \(E\) fields.

    3. Compute the linear transfer matrices (\(M_x\) and \(M_y\)) through a bend-only region of length \(\ell_b = 51/360\cdot 2\pi R_0\).

    4. Compute the linear transfer matrices (\(M_x\) and \(M_y\)) through a bend-and-focus region of length \(\ell_q = 39/360\cdot 2\pi R_0\).

    5. Compute the two products (once for \(x\) and once for \(y\)) of the relevant matrices above, representing horizontal and vertical motion through 1/4 of the ring and also compute the corresponding matrices for one complete revolution.

    6. Using the last two matrices above, estimate the horizontal and vertical betatron tunes of the complete ring. Additionally, estimate the values of \(\beta_x\) and \(\beta_y\) from these matrices.

# a
c = 2.99792458e8
R0 = 7.112                              # m
B0 = (pmag[2]/1000)*10/c*1e8/R0    # T
Brho = B0*R0                            # T-m
B0 # Tesla
## [1] 1.451295
# b (i)
Kxb = 1/R0^2
Kyb = 0
print(paste("Kxb =",round(Kxb,4),"/m^2;    Kyb =",
   round(Kyb,4),"/m^2"))
## [1] "Kxb = 0.0198 /m^2;    Kyb = 0 /m^2"
# b (ii)
Vq = 20e3                               # V
aq = 50e-3                              # m
dEdx = 2*Vq/aq^2                        # V/m^2
Kxq = 1/R0^2 - dEdx/c/Brho
Kyq =          dEdx/c/Brho
print(paste("Kxq =",round(Kxq,4),"/m^2;    Kyq =",
   round(Kyq,4),"/m^2")) 
## [1] "Kxq = 0.0146 /m^2;    Kyq = 0.0052 /m^2"
# c
sb = 51/360*2*pi*R0
MxB = matrix( c(
      cos(sqrt(Kxb)*sb),           1/sqrt(Kxb)*sin(sqrt(Kxb)*sb),
  -sqrt(Kxb)*sin(sqrt(Kxb)*sb),           cos(sqrt(Kxb)*sb)    ), ncol=2, byrow=TRUE)
MyB = matrix( c(
             1,   sb,
             0,    1      ), ncol=2,byrow=TRUE)
MxB
0.6293204 5.5270621
-0.1092725 0.6293204
MyB
1 6.330519
0 1.000000
# d 
sq = 39/360*2*pi*R0
MxQ = matrix( c(
      cos(sqrt(Kxq)*sq),           1/sqrt(Kxq)*sin(sqrt(Kxq)*sq),
  -sqrt(Kxq)*sin(sqrt(Kxq)*sq),           cos(sqrt(Kxq)*sq)     ), ncol=2, byrow=TRUE)
MyQ = matrix( c(
      cos(sqrt(Kyq)*sq),           1/sqrt(Kyq)*sin(sqrt(Kyq)*sq),
  -sqrt(Kyq)*sin(sqrt(Kyq)*sq),           cos(sqrt(Kyq)*sq)       ),   ncol=2,byrow=TRUE)

MxQ
0.8337492 4.5696149
-0.0667151 0.8337492
MyQ
0.9400210 4.743806
-0.0245289 0.940021
# e
Mx4 = MxQ %*% MxB
My4 = MyQ %*% MyB

Mx4
0.0253622 7.483936
-0.1330910 0.155957
My4
0.9400210 10.69463
-0.0245289 0.78474
Mx = Mx4 %*% Mx4 %*% Mx4 %*% Mx4
My = My4 %*% My4 %*% My4 %*% My4

Mx
0.9580772 -2.6693490
0.0474705 0.9114971
My
-0.3943453 17.9808512
-0.0412405 -0.6554188
# f
nux0 = acos(sum(diag(Mx))/2)/2/pi
nuy0 = acos(sum(diag(My))/2)/2/pi

nux0
## [1] 0.05779504
nuy0
## [1] 0.3379453
# May be better to take tune of 1/4 of ring, multiply by 4:
nux = 4*acos(sum(diag(Mx4))/2)/2/pi
nuy = 4*acos(sum(diag(My4))/2)/2/pi

nux
## [1] 0.942205
nuy
## [1] 0.3379453
betax = Mx[1,2]/sin(2*pi*nux)  #  m
betay = My[1,2]/sin(2*pi*nuy)  #  m

betax
## [1] 7.514882
betay
## [1] 21.12474
  1. Admittance/Emittance. Assume a round aperture of radius \(a_0\) = 45 mm, and constant values of \(\beta_x\) = 7.5 m, \(\beta_y\) = 21 m and \(D_x\) = 8 m for the Muon Storage Ring.

    1. Estimate the horizontal and vertical transverse admittances \(A_{x0}\), \(A_{y0}\) for the central momentum (\(\delta \equiv \delta p/p\) =0); i.e., what is the maximum phase space area for each dimension in which the beam can reside?

    2. If the incoming beam is Gaussian with rms emittance 20 \(\pi\) mm-mr in each degree of freedom (\(x\) and \(y\)), estimate the maximum fraction that can be contained within the overall admittance for the central momentum?

    3. In reality, the incoming beam is squeezed horizontally to fit through the inflector, which has a horizontal full-width aperture of 18 mm. How does this affect (quantitatively) your answer in (b) above?

    4. The incoming beam has a momentum spread of \(\delta_{rms}\) = 2% (Gaussian). What is the maximum momentum acceptance (\(\hat{\Delta p}/p\)) of the ring? Estimate what fraction of the incoming distribution might possibly circulate in the ring.

    5. Estimate the overall fraction of incoming particles that will survive long-term (i.e., many turns after injection, but neglecting decays) in the Storage Ring.

    6. What other factors might affect the overall storage efficiency?

a0 = 45
betx = 7.5
bety = 21
Dx = 8
eps_rms = 20

# a
Ax0 = a0^2/betx
Ay0 = a0^2/bety
Ax0
## [1] 270
Ay0
## [1] 96.42857
# b
f_storX = ifelse(Ax0/6/eps_rms >1,
                  1, Ax0/6/eps_rms)
f_storY = ifelse(Ay0/6/eps_rms >1,
                  1, Ay0/6/eps_rms)
f_stor = f_storX * f_storY
f_stor
## [1] 0.8035714
# c) squeezed injection
#       pi betx*6*eps_rms --> 18*dxInj
dxInj = pi*betx*6*eps_rms/18
dxInj
## [1] 157.0796
f_sqzX = (2*a0)/dxInj
f_sqz  = f_sqzX*f_storY

# d
del_rms = 0.02
dpMax = a0/(Dx*1000)
dpMax
## [1] 0.005625
f_mom = 2*dpMax/sqrt(2*pi)/del_rms
f_mom
## [1] 0.224405
# e
f_sqz*f_mom
## [1] 0.1033189
# f) Admittance vs. dp/p goes to zero:  another factor of two, say?
f_sqz*f_mom/2
## [1] 0.05165944

Tuesday (due Wednesday)

  1. Pion Decay. A beam of positive pions (\(\pi^+\); rest mass 140 MeV/c\(^2\)) is accelerated, so it has momentum \(p_{\pi^+}\) = 3.1 GeV/c in the laboratory frame. A \(\pi^+\) particle decays into a positive muon (rest mass 105 MeV/c\(^2\)) and a muon neutrino (\(m_{\nu_\mu}\) = 0 MeV/c\(^2\)).

    1. What is the energy of the muon in the rest frame of the \(\pi^+\) particle?

    2. What is the energy of the muon in the laboratory frame if it is found traveling in the beam direction?

    3. What is the energy of the muon in the laboratory frame if it is found traveling at 0.25 rad with respect to the beam direction?

# a) In the pion rest frame:
mc2pi = 139.57 # MeV
mc2mu = 105.65 # MeV
Estar = (mc2pi^2+mc2mu^2)/2/mc2pi
Estar
## [1] 109.7718
pstar = (mc2pi^2 - mc2mu^2)/2/mc2pi
pstar
## [1] 29.79817
# b) In the lab frame:
#    E = gamma_pi * [Estar + beta_pi * pstar * cos(thetastar)]
#       cos(thetastar) = cos(0) = 1
gamma_pi = sqrt(3100^2+mc2pi^2)/mc2pi
gamma_pi
## [1] 22.23358
beta_pi  = sqrt(1-1/gamma_pi^2)
beta_pi
## [1] 0.998988
E = gamma_pi*(Estar + beta_pi * pstar * cos(   0))
E/1000 # GeV
## [1] 3.10247
# c)  If thetastar = 0.25 rad, ...
E = gamma_pi*(Estar + beta_pi * pstar * cos(0.25))
E      # GeV
## [1] 3081.895
  1. Emittance Measurement. A plane of horizontal wires is used to detect charged particles passing through it to create a rough measurement of the vertical particle distribution at that location. A beam of muons yields a roughly Gaussian vertical distribution with rms \(\sigma_1\) = 12 mm. Located downstream is a second set of wires. This set is translated longitudinally until the rms beam size at its location is minimized. Under this condition, \(\sigma_2\) = 3.2 mm. If the two sets of wires are located 1.4 m apart when these two measurements are made, and if there are no gradient fields in the region between them (and no dispersion),

    1. What is the vertical transverse emittance of the muon beam?

    2. What are the values of the Courant-Snyder parameters \(\alpha\) and \(\beta\) at the two wire locations?

Solution

Use \(x_2 = x_1 + x'\), where \(x' = x'_1=x'_2 = (x_2-x_1)/L\), and also we know that \(\langle x_2x'_2\rangle = 0\) since at a waist at location 2. So we concentrate on location 2, and \[ \langle x_2'^2\rangle = \langle x'^2\rangle = \langle (\frac{x_2-x_1}{L})^2\rangle = \frac{\sigma_2^2 + \sigma_1^2 - 2\langle x_1x_2\rangle}{L^2}, ~~~{\rm and} \\ \langle x_2x'_2\rangle = \langle \frac{x_2(x_2-x_1)}{L} \rangle = \frac{\sigma_2^2-\langle x_1x_2\rangle}{L} = 0 \] from which \(\langle x_1x_2\rangle = \sigma_2^2\). Thus, \[ \langle x_2'^2\rangle = \langle x'^2\rangle = \frac{\sigma_2^2 + \sigma_1^2 - 2\sigma_2^2}{L^2} = \frac{\sigma_1^2 - \sigma_2^2}{L^2} \] So, \[ \epsilon = \pi\sqrt{\sigma_2^2\cdot\frac{\sigma_1^2-\sigma_2^2}{L^2} - 0} = \pi\frac{\sigma_2}{L}\sqrt{\sigma_1^2-\sigma_2^2} \]

L    =  1.4 # m   
sig1 = 12   # mm
sig2 =  3.2 # mm
eps = sig2/L*sqrt(sig1^2-sig2^2)
bet2 = sig2^2/eps
alf2 = 0
bet1 = sig1^2/eps
alf1 = (sig1^2-sig2^2)/L/eps
out = cbind.data.frame(c(bet1,bet2),c(alf1,alf2),row.names=c("Loc1","Loc2"))
print(paste("emittance =",round(eps,2),"pi mm-mr"))
## [1] "emittance = 26.44 pi mm-mr"
kable(round(out,3),col.names = c("beta","alpha"),row.names=TRUE)
beta alpha
Loc1 5.447 3.614
Loc2 0.387 0.000
  1. Time of Flight. Suppose an ideal particle takes time \(T_0\) to travel an ideal distance \(L_0\) at an ideal speed \(v_0\). A particle with speed \(v = v_0(1+dv/v_0)\) will take a different amount of time \(T\) for it to reach the same longitudinal location, both due to its speed and due to the fact that it may take a different path to get there.

    1. Show that the relative differential time will be \(dT/T_0 = dL/L_0 - dv/v_0\).

    2. Show that, in terms of momentum, \(dv/v_0 = (1/\gamma_0^2) dp/p_0\).

    3. When traveling through bending magnets, particles will take different paths, determined by the dispersion function \(D(s)\). Show that \[ dL = \frac{dp}{p_0} \int \frac{D(s)}{\rho(s)} ds \] where \(\rho\) is the local bending radius.

    4. For a storage ring, show that \[ dT/T_0 = ( \alpha_p - 1/\gamma_0^2 ) ~dp/p_0 \] where \(\alpha_p = \langle D/\rho \rangle\) where the average is taken over the ring circumference. Note that \(\rho\) is infinite outside of the dipole magnets

    5. For a muon beam with a total momentum spread (\(5\sigma\)) of 7.5% at a central momentum of 3.09 GeV/c, and a ring with \(\langle D \rangle\) = 1.5 m within the dipoles magnets, \(\rho\) = 17.5 m, and a packing fraction (percent of circumference with bending) of 22%, determine the spread in revolution times for the distribution.

    6. If the particles all arrived simultaneously to the ring, what will be the total time spread of the distribution after 4 revolutions?

    7. How does this time spread compare to the bunch length of the beam coming from the Recycler?

# e)
Dave = 1.5
rho  = 17.5
Brho = 10/2.9979*3.09
pfrc = 0.22
gammu = 3.09/0.10566

eta = Dave/rho*pfrc - 1/gammu^2
dT = (2*pi*rho/pfrc/2.9979e8)*eta*0.075 
dTtot = dT*1e9   # ns  (total spread)
dTtot
## [1] 2.211646
# f) Now, for 4 turns:
4*dTtot
## [1] 8.846585
# g) The g-2 ring revolution period:
2*pi*7.112/2.9979e8* 1e9  # ns
## [1] 149.0577
  1. Quad Multipoles. The electric field can be expanded in Cartesian coordinates in terms of multipoles \(b_n\) and \(a_n\) according to \[ E_x-iE_y = (b_n-ia_n)\frac{(x+iy)^n}{r_0^n} \] where \(E_z = 0\).

    1. Show that the equation above satisfies Maxwell’s equations order by order.
    2. Show that, in cylindrical coordinates where \(r = x+\rho_0\) and \(z=y\), the expression for the electric field violates Maxwell at lowest order \(n=1\).
    3. We can however write \(E_x\) in the midplane (\(z=0\)) in the form \[ E_x = \sum_n b_n\frac{x^n}{r_0^n} \] with \(x=\rho-\rho_0\), and \(b_n\) as given in the table below for \(r_0\) = 0.045 m and plate voltage of 27.2 kV. Compute the contribution to the amplitude dependent tune shift from each of the multipoles. Compare the total tune shift with the result from tracking.
Multipole (n) Value [\({\rm V/m}\)]
1 \(1.008\times 10^{6}\)
2 \(-7.043\times 10^{4}\)
3 \(-7.048\times 10^{3}\)
4 \(-8.561\times 10^{2}\)
5 \(2.520\times 10^{3}\)
6 \(2.605\times 10^{2}\)
7 \(-1.444\times 10^{3}\)
8 \(-1.830\times 10^{2}\)
9 \(-9.852\times 10^{4}\)
10 \(7.400\times 10^{3}\)
11 \(1.189\times 10^{3}\)
12 \(-4.795\times 10^{1}\)
13 \(2.552\times 10^{3}\)

Wednesday (due Thursday)

  1. Finish morning Computer Session’s assignment for Thursday.

  2. Fast Rotation. Imagine that the initial distribution has zero emittance, zero energy spread and zero bunch length. The particles share a common revolution period \(T\) and the time dependence of the intensity signal at a fixed point in the ring (a fiber harp for example), is \[ I(t) = \sum_{n=0}^\infty\delta(t-nT) \]

    1. Compute the Fourier transform of the intensity to determine the frequency spectrum.

    2. If a particle has a momentum offset \(\Delta\), the revolution period \(T\rightarrow T(1+\Delta)\) and \[ I(t,\Delta) = \sum_{n=0}^\infty\delta(t-nT(1+\Delta)). \] So, if \(\rho(\Delta)\) is the normalized momentum distribution then the signal in the detector is \[ S(t) = \sum_{n=0}^\infty\int\rho(\Delta)\delta(t-nT(1+\Delta))d\Delta. \] Suppose the momentum distribution is Gaussian with width \(\Delta_0\).

      1. Integrate the equation for \(S(t)\) to find an expression for \(S(t)\) in terms of \(\Delta_0\) and \(T\). (This is the Fast Rotation signal that we measure.)

      2. Compute the Fourier transform \(F(\omega)\) of \(S(t)\).

      3. Plot \(F(\omega)\) over the frequency range defined by the collimators. How is the width of the frequency distribution related to the \(\Delta_0\)?

      4. In the experiment, our measurement of \(S(t)\) starts at some time \(t=t_0\) that excludes the first few microseconds (turns) of the fill. How is the frequency distribution affected?

  3. Messing with Perfection. Imagine a perfectly constructed Muon \(g\)-2 Storage Ring.

    1. Treat one of the “long” quadrupoles as a thin lens. What would be its focal length for a quad voltage of 18.3 kV? Is the thin lens treatment valid?

    2. If one of the long quadrupoles is displaced (as a unit) vertically upward by 2 mm, make a plot of the resulting vertical closed orbit distortion.

    3. If the same quad is aligned perfectly vertically, but displaced radially outward by 2 mm, plot the horizontal closed orbit distortion that results.

    4. Using Bmad, verify your results for parts b) and c).

    5. Using Bmad, give all of the quadrupole “units” (long and short) random displacements, horizontal and vertical, chosen from a Gaussian distribution with mean of zero and rms of 2 mm. Plot the overall closed orbit distortions (horizontal and vertical) created.

# a)
R0 = 7.112                              # m
sq = 39/360*2*pi*R0*2/3
c = 2.99792458e8
Vq = 18.3e3                             # V
aq = 50e-3                              # m
B0 = (3094/1000)*10/2.99792458/R0    # T
Brho = B0*R0                            # T-m
dEdx = 2*Vq/aq^2                        # V/m^2
K = dEdx/c/Brho
Finv = K*sq
F = 1/Finv
F  # m
## [1] 65.48424
# if length of the quad << F, then "thin"
sq/F
## [1] 0.04928397
# b)
thetay = 2/F # mr
psiy = c(0:24)*2*pi*nuy/24
dy = betay*thetay/2/sin(pi*nuy)*cos(psiy-pi*nuy)
plot(psiy/2/pi/nuy*2*pi*R0,dy,ylim=c(-1,1),typ="l")
abline(h=0)

# c)
thetax = -2/F
psix = c(0:24)*2*pi*nux/24
dx = betax*thetax/2/sin(pi*nux)*cos(psix-pi*nux)
plot(psix/2/pi/nux*2*pi*R0,dx,ylim=c(-1,1),typ="l")
abline(h=0)

  1. Full Momentum. Suppose that muons were injected into the Muon g-2 Storage Ring such that (a) the entire acceptance was filled uniformly (in \(x\), \(x'\), \(y\), \(y'\)) and that the momentum of each particle is also random with uniform distribution between \(\pm 6\times 10^{-3}\) of the magic momentum. Remember that the Storage Ring has a circular aperture of radius \(r_0\) = 45 mm.

    1. Compute the final momentum distribution of particles that can survive long-term.

    2. What is the rms of this final momentum distribution?

a0 = 45   # mm
D  = 8100 # mm
Npar = 10000

phix0 =     runif(Npar,0,1)*2*pi
phiy0 =     runif(Npar,0,1)*2*pi
ax0  = sqrt(runif(Npar,0,1))*a0
ay0 = sqrt(runif(Npar,0,1))*a0
x0  = ax0*cos(phix0)
px0 = ax0*sin(phix0)
y0  = ay0*cos(phiy0)
py0 = ay0*sin(phiy0)
del = runif(Npar,-1,1)*0.006
df = data.frame(cbind(x0,px0,y0,py0,del))
df2 = subset(df,x0^2+y0^2<a0^2)
length(df2$x0)
## [1] 9112
plot(df2$x0,df2$px0,asp=1,pch=".")

plot(df2$y0,df2$py0,asp=1,pch=".")

plot(df2$x0,df2$y0,asp=1,pch=".")

hist(df2$del,breaks=25)

df3 = subset(df2, y0^2+py0^2 + (D*abs(del)+sqrt((x0-D*del)^2+px0^2) )^2 < a0^2)
length(df3$x0)
## [1] 1259
plot(df3$x0,df3$px0,asp=1,pch=".")

plot(df3$y0,df3$py0,asp=1,pch=".")

plot(df3$x0,df3$y0,asp=1,pch=".")

#  a)
hist(df3$del,breaks = 25)

#  b)
print(paste( "rms dp/p =",round(sd(df3$del)*100,3),"%"))
## [1] "rms dp/p = 0.151 %"
  1. Dispersion Matching. The M5 beam line bringing muons to the Muon g-2 Storage Ring has both horizontal and vertical dispersions of zero, by design, in order to “squeeze” the beam through the inflector at the injection point. We saw earlier that this results in a relatively low injection efficiency. Suppose that the beam delivery system is redesigned to match the dispersion of the ring upon injection. And suppose that the inflector aperture could be redesigned to accommodate this change. Assume that the incoming beam has rms emittances \(\epsilon_x\) = \(\epsilon_y\) = 15 \(\pi\) mm-mr, and rms momemtum spread of \(12\times 10^{-3}\) (\(\Delta p/p\)) (Gaussian distributions for all variables). (Also, assume that the kicker system works as it should.)

    1. Estimate the fraction of the incoming muons that would you expect to capture upon injection using this new system.
    2. Estimate the resulting rms momentum spread of the particles that can survive long-term.
    3. What would be the result if the aperture were changed into a 45 mm \(\times\) 45 mm square, rather than a circle?
ap   = 45   # mm (aperture)
Dx   = 8100 # mm (dispersion)
eps  = 15   # pi mm-mr
sigp = 0.012
betax = 7.6
betay =  21
# Emittance in normalized phase space
maxAdm = 6*(pi*eps)*betax
# Admittance of particle with momentum del=x:
Adel = function(x){ pi*(ap-Dx*abs(x))^2 }
curve(Adel,-ap/Dx,ap/Dx)
abline(h=maxAdm,col="blue")

Npar = 500000
del = rnorm(Npar,0,sigp)
x   = rnorm(Npar,0,sqrt(betax*eps)) + Dx*del
PX  = rnorm(Npar,0,sqrt(betax*eps))
y   = rnorm(Npar,0,sqrt(betay*eps))
PY  = rnorm(Npar,0,sqrt(betay*eps))
a0 = sqrt( (x-Dx*del)^2 + PX^2 )
b0 = sqrt(  y^2         + PY^2 )
plot(x,PX,pch=".",asp=1)

df = data.table(cbind(x,PX,y,PY,del,a0,b0))
capt0 = df[Dx*abs(del) < ap]
capt = capt0[b0   <   sqrt( ap^2 - (Dx*del)^2 )  & 
                     b0^2 + (Dx*abs(del) + a0)^2 < ap^2  ]
plot(capt$x,capt$PX,pch=".",asp=1)

hist(capt$del)

# a)
capt[,.N]/Npar
## [1] 0.19117
# b)
sd(capt$del)
## [1] 0.002067055
# c)  change to a square...
captsq = capt0[b0 < ap   &   Dx*abs(del)+a0 < ap  ]
plot(captsq$x,captsq$PX,pch=".",asp=1)

hist(captsq$del)

captsq[,.N]/Npar
## [1] 0.245308
sd(captsq$del)
## [1] 0.002376199

Thursday (Group Presentations due Friday)

  1. Finish morning Computer Session’s assignment for Friday.
  2. FINAL topic 1:
  3. FINAL topic 2:
  4. FINAL topic 3: