Homework 6 Discussion

Some General Comments:

• grade: 30 possible
  + 10 pts each problem (+10 for bonus)
  + ave score: 27/30 (incl. bonus)

• use of data.table, etc.
  • df2 = df[PDGid == -13]; nrow(df2)
    

• some leeway given for Prob. 1c (vague wording)

Muon Decay Kinematics

A beam of \(\pi^+\) particles (\(m_{\pi}\) = 140 MeV/\(c^2\)) is accelerated to a momentum of \(p_{\pi_+}\) = 3.1 GeV/\(c\) in the laboratory frame. A \(\pi^+\) particle decays into a positive muon (\(m_{\mu}\) = 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^+\)?
   
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 and angle of 250 mrad 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, ...
#  Note:  wording was a bit ambiguous; graded generously
E = gamma_pi*(Estar + beta_pi * pstar * cos(0.25))
E/1000 # GeV

[1] 3.081895

G4beamline Problem (optional)

1. Assume that at the downstream end of the lithium lens (\(z\) = 0.577 m) we have a collection of pions all at 3.094 GeV/\(c\). Based on the exponential decay law, how many pions will survive at the end of M3 (\(z\) = 279.230 m) in the lab frame?

mc2Pi = 139.57  # MeV
Epi   =  3094   # MeV
tauPi = 2.6e-8  # s
gamPi = Epi/mc2pi
tauLab = gamPi*tauPi # s
dt2M3 = (279.23 - 0.577)/2.9979e8  # s
Nf = exp(-dt2M3/gamPi/tauPi)
Nf

[1] 0.199355

The pion lifetime in its own frame is \(\tau\) = \(2.6\times 10^{-8}\) s. In the lab frame the pion lifetime is \(\gamma\tau\) = \(5.7637028\times 10^{-7}\) where \(\gamma\) = 22.1680877 for 3.094 GeV/c pions. The time needed to reach the end of M3 is dt2M3 = \(9.2949398\times 10^{-7}\) s. The exponential decay law predicts about 20% of the pions remain at the of M3, or 80% of the pions have decayed.

2. Now let’s open the G4beamline GUI and download and open the G4beamline file HW1_G4_M2M3.g4bl. This is a beam line description of the “M2-M3” beam line at the Fermilab Muon Campus; the beam line starts just in front of a target used for producing muon beams. Start to play with the visualization tools. Note: A demonstration of using G4beamline was provided in class.

[Simply run G4beamline and proceed to next question…]

3. Starting at the downstream end of the lithium lens, the G4beamline code (HW1_G4_M2M3.g4bl) tracks a distribution of 3000 pions (param nparticles = 3000) with a momentum of 3.1 GeV/\(c\). Running the code should take about 10 minutes or less, depending on the speed of your computer. The distribution at the end of M3 should be generated by the code, and written to a text file called HW1_particles_EndM3.txt. Open the file and exam with an appropriate editor and/or program (such as R). How many pions will survive to the end of M3? How does it compare with the exponential decay law? Note: PDGid for pions is 211 and PDGid for muons is -13.

Running the simulation predicts 573 pions at the end of M3. Starting with 3000 pions, the exponential decay calculation above would predict \(3000~\times\) Nf = 598, which is very close.

4. For the distribution you have created for Part c above, make a histogram of the muon and pion momentum distributions. What do you observe?

In addition to pions we have \(\approx 100\) muons. The pion and muon momentum histograms are created and shown below:

library(data.table)

data.table 1.12.2 using 2 threads (see ?getDTthreads).  Latest news: r-datatable.com

df = fread("HW1_particles_EndM3.txt",skip = 2)
head(df)

colnames(df[,8]) = "P0"
dfpi = df[V8 == 211]
dfmu = df[V8 == -13]
nrow(dfpi)  # number of pions

[1] 610

nrow(dfmu)  # number of muons

[1] 84

par(mfrow=c(1,2))
hist(dfpi$V6,breaks=50,main="pion distribution",xlab="pion energy [MeV]")
hist(dfmu$V6,breaks=50,main="muon distribution",xlab="muon energy [MeV]")

Beam Debunching

Protons transported to the Fermilab Booster ring are first accelerated through a linear accelerator. The particles leaving this “linac” are arranged in bunches spaced 5 ns apart, where each bunch has an average kinetic energy of 400 MeV, a relative momentum spread of 0.3% and a time spread of 0.5 ns (rms values for each; assume Gaussian distributions). After leaving the linac, the particles travel approximately 100 m and then go through a cavity. The cavity is operating sinusoidally with a 200 MHz frequency and can be used to lower the overall momentum spread of the beam. The cavity is located approximately 100 m before the Booster injection point.

1. Write a short code that transports several thousand particles within one bunch, with the initial conditions given above, through 100 m of drift, then through the cavity, then 100 m further. Plot the longitudinal phase space distribution (\(\Delta p/p\) vs. \(\Delta t\)) at the beginning and end of the transport system, as well as just before and after the cavity. Treat the cavity as a “thin” device for our purposes, but retain its sinusoidal time structure.

2. What overall effective voltage is required of the cavity to minimize the momentum spread of the beam? What minimum value do you arrive at?

Np    = 10000
Veff  = 2.5e6
gamma = (938+400)/938
beta  = sqrt(1-1/gamma^2)
dt0   = rnorm(Np,0,0.05)    # ns
del   = rnorm(Np,0,0.002)  # dp/p
par(mfrow=c(2,2))
plot(dt0,del*100,xlim=c(-3,3),ylim=c(-1,1)*1.2,pch=".",
   xlab="dt [ns]",ylab="dp/p [%]",main="exit of linac")
text(1.7,0.8,paste("dp/p rms = ",round(sd(del)*100,2)," %" ))

dt1   = dt0 -(100/2.9979e8*1e9)/beta/gamma^2*del
plot(dt1,del*100,xlim=c(-3,3),ylim=c(-1,1)*1.2,pch=".",
   xlab="dt [ns]",ylab="dp/p [%]",main="entrance to cavity")
curve(sin(2*pi*200e6*x/1e9),add=TRUE,col="blue",lty=2)

del2  = del + sin(2*pi*200e6*dt1/1e9)*Veff/400e6*gamma/(1+gamma)
plot(dt1,del2*100,xlim=c(-3,3),ylim=c(-1,1)*1.2,pch=".",
   xlab="dt [ns]",ylab="dp/p [%]",main="exit of cavity")
text(1,0.8,paste("Veff =", round(Veff/1e6,2),"MV\n","rms = ",round(sd(del2)*100,3)," %" ))

dt2   = dt1 -(100/2.9979e8*1e9)/beta/gamma^2*del2
plot(dt2,del2*100,xlim=c(-3,3),ylim=c(-1,1)*1.2,pch=".",
   xlab="dt [ns]",ylab="dp/p [%]",main="entrance to Booster")
text(1,0.8,paste("Veff =", round(Veff/1e6,2),"MV\n","rms = ",round(sd(del2)*100,3)," %" ))

Adiabatic Ramping

Synchrotron motion in the longitudinal degree of freedom can be characterized by the following set of difference equations:

SynchTrk = function(U){
   U[2] = U[2] + k*sin( 2*pi*(U[1]-us) )
   U[1] = U[1] - g*U[2]
   return(U)
}

where u and v are “energy” and “phase” variables, and k and g are constants determined by the accelerating voltage and time of flight between accelerating stations. The synchronous phase is given by the parameter us. The above function SynchTrk will hence track the u-v variables through one turn about a synchrotron or, say, between two successive cavities in a linac.

g  = 0.02
k  = 0.005
us = 0

u = 0.01
v = 0

1. Create a small code chunk that uses the above function to track a single particle through nturn iterations. From your code, determine the approximate period of a synchrotron oscillation for values of k = 0.005, g = 0.02 and us = 0.0833333. [A good initial condition for your particle is (u,v) = (0.01, 0).]

Tip: In R, employ the following plot attributes to make consistent looking plots:

# prior to the loop:
plot(U[1],U[2],xlim=c(-1,1)*0.15,ylim=c(-1,1)*0.05,asp=1,typ="n")
#     asp is the "aspect ratio"; typ="n" means do not plot anything yet

#... within the loop:
   points(U[1],U[2],pch=".")    # adds points to the previous plot
#...

nturn =  5000
U = c(u,v)
plot(u,v,xlim=c(-1,1)*0.15,ylim=c(-1,1)*0.05,asp=1,typ="n")

i = 0
while(i<250){  # the phase space trajectory appears to "close" after about 250 turns
   i = i+1
   U = SynchTrk(U)
   points(U[1],U[2],pch=".")
}

2. Next, we wish to adjust the synchronous phase from an initial value of us0 = 0 (zero degrees) to a value of usmax = \(\pi/6\) (30 degrees). In order to explore various “ramp rates” of this phase, create a function usR(n) that
   1. maintains usR = us0 until n = on,
      
   2. varies usR linearly from us0 to usmax between turns n = on and n = off,
      
   3. maintains usR = usmax for n > off.
      Make a plot of your function over the range 0 < n < 400 turns using values of on = 100, off = 300.

on    =   100
off   =   300
us0   =     0
usmax =   1/12

usR = function(x){
   ifelse(x<on,us0,ifelse(x>=off,usmax,
      (usmax-us0)/(off-on)*(x-on) ))
}
curve(usR,0,400)

3. Modify your tracking code to employ your phase ramp function in the calculation. Make an intial test by setting usmax = 0 and on = 1, off = 2 then re-tracking for a sufficient number of turns to make sure you get the same result as previously.

nturn =  400
U = c(u,v)
plot(u,v,xlim=c(-1,1)*0.15,ylim=c(-1,1)*0.05,asp=1,typ="n")

i = 0
usmax = 0
on = 1
off = 2
while(i<nturn){  # the phase space trajectory appears to "close" after about 250 turns
   i = i+1
   us = usR(i)
   U = SynchTrk(U)
   points(U[1],U[2],pch=".")
}

nturn =  8000
us0   =     0
usmax =   1/12
on    =   250
offJ  = c(350, 400, 600, 2000, 3000, 6000)

4. Finally, set usmax = 0.0833333, on = 250 and make phase space plots for various “ramp off” times: off = 350, 400, 600, 2000, 3000, 6000. For consistency, use nturn = 8000 for each case. Comment on the results.

par(mfrow=c(2,3))
j = 0
while(j<6){
   j = j+1
   off = offJ[j]
   U = c(u,v)
   plot(u,v,xlim=c(-1,1)*0.15,ylim=c(-1,1)*0.05,asp=1,
      typ="n",main=paste("off =", off))
   i = 0
   while(i<nturn){
      i = i+1
      us = usR(i)
      U = SynchTrk(U)
      points(U[1],U[2],pch=".",   # add colors to envision times during ramping
         col=ifelse(i<on,"green",ifelse(i>off,"red","black") ))
   }
   abline(v=c(0,usmax),lty=3,col="brown")
}

---

1. Northern Illinois University and Fermi National Accelerator Laboratory↩