diff --git a/simc/shared/musc.f b/simc/shared/musc.f
index 37d680d3c9020e20cfc2c1e9840831d8085ad1bc..2fd7dac6abedb2d7c982d8dd32c09eb2e9d26d43 100644
--- a/simc/shared/musc.f
+++ b/simc/shared/musc.f
@@ -25,34 +25,21 @@ C-_____________________________________________________________________
implicit none
- real*8 Es, epsilon
- parameter (Es = 13.6) !MeV
- parameter (epsilon = 0.088)
+ real*8 musc_pdg, musc_with_tail
real*8 rad_len, dth, dph
- real*8 beta, theta_sigma
- real*8 m2, p
-
- real*8 nsig_max
- parameter(nsig_max=99.0e0) !max #/sigma for gaussian ran #s.
-
- real*8 gauss1
-
-! Compute scattering angles, THETA_SCAT from a gaussian distribution,
-! PHI_SCAT from uniform distribution.
+ real*8 beta, m2, p
if (rad_len.eq.0) return
if (p.lt.25.) write(6,*)
> 'Momentum passed to musc.f should be in MeV, but p=',p
+! Compute new trajectory angles (units are rad)
+
beta = p / sqrt(m2+p*p)
-c theta_sigma = Es/p/beta * sqrt(rad_len) * (1+epsilon*log10(rad_len))
-c Better form for beta .ne. 1
- theta_sigma = Es/p/beta * sqrt(rad_len) * (1+epsilon*log10(rad_len/beta**2))
-! Compute new trajectory angles (units are rad)
+ dth = dth + musc_with_tail(beta, p, rad_len)
+ dph = dph + musc_with_tail(beta, p, rad_len)
- dth = dth + theta_sigma * gauss1(nsig_max)
- dph = dph + theta_sigma * gauss1(nsig_max)
return
end
diff --git a/simc/shared/musc_ext.f b/simc/shared/musc_ext.f
index ddc78414df7843bef4c519d282e837d3b23559ff..d44b108bc9db2409dd59b72d9c4bd844ad28ba8c 100644
--- a/simc/shared/musc_ext.f
+++ b/simc/shared/musc_ext.f
@@ -15,13 +15,9 @@ C-_____________________________________________________________________
parameter (epsilon = 0.088)
real*8 rad_len, x_len, dth, dph, x, y
- real*8 beta, g1, g2, theta_sigma
+ real*8 beta, rx1, rx2, ry1, ry2
real*8 m2, p
-
- real*8 nsig_max
- parameter(nsig_max=99.0e0) !max #/sigma for gaussian ran #s.
-
- real*8 gauss1
+ real*8 musc_pdg, musc_with_tail
if (rad_len.eq.0) return
if (x_len.le.0 .or. rad_len.lt.0) then
@@ -33,22 +29,18 @@ C-_____________________________________________________________________
if (p.lt.10.) write(6,*)
> 'Momentum passed to musc_ext.f should be in MeV, but p=',p
- beta = p / sqrt(m2+p*p)
-c theta_sigma = Es/p/beta * sqrt(rad_len) * (1+epsilon*log10(rad_len))
-C Better form for beta .ne. 1
- theta_sigma = Es/p/beta * sqrt(rad_len) * (1+epsilon*log10(rad_len/beta**2))
-
! Compute new trajectory angles and displacements (units are rad and cm)
+ beta = p / sqrt(m2+p*p)
- g1 = gauss1(nsig_max) ! gaussian, truncated at 99 sigma
- g2 = gauss1(nsig_max)
- dth = dth + theta_sigma*g1
- x = x + theta_sigma*x_len*g2/sqrt(12.) + theta_sigma*x_len*g1/2.
+ rx1 = musc_with_tail(beta, p, rad_len)
+ rx2 = musc_with_tail(beta, p, rad_len)
+ ry1 = musc_with_tail(beta, p, rad_len)
+ ry2 = musc_with_tail(beta, p, rad_len)
- g1 = gauss1(nsig_max) ! gaussian, truncated at 99 sigma
- g2 = gauss1(nsig_max)
- dph = dph + theta_sigma*g1
- y = y + theta_sigma*x_len*g2/sqrt(12.) + theta_sigma*x_len*g1/2.
+ dth = dth + rx1
+ x = x + rx2*x_len/sqrt(12.) + rx1*x_len/2.
+ dph = dph + ry1
+ y = y + ry2*x_len/sqrt(12.) + ry1*x_len/2.
return
end
diff --git a/simc/shared/musc_pdg.f b/simc/shared/musc_pdg.f
new file mode 100644
index 0000000000000000000000000000000000000000..aba861b7f34470ea9f6b2be5e344e45af9b6db3c
--- /dev/null
+++ b/simc/shared/musc_pdg.f
@@ -0,0 +1,28 @@
+ real*8 function musc_pdg(beta,p,rad_len)
+C+_____________________________________________________________________
+!
+! MUSC - Implementation of the PDG formula for multiple scattering
+!
+! According to Particle Data Booklet, July 1994
+!
+C-_____________________________________________________________________
+
+ implicit none
+
+ real*8 Es, epsilon
+ parameter (Es = 13.6) !MeV
+ parameter (epsilon = 0.088)
+
+ real*8 rad_len
+ real*8 beta, theta_sigma, p
+
+ real*8 nsig_max
+ parameter(nsig_max=99.0e0) !max #/sigma for gaussian ran #s.
+
+ real*8 gauss1
+
+ theta_sigma = Es/p/beta * sqrt(rad_len) * (1+epsilon*log10(rad_len/beta**2))
+
+ musc_pdg = gauss1(nsig_max) * theta_sigma
+ return
+ end
diff --git a/simc/shared/musc_with_tail.f b/simc/shared/musc_with_tail.f
new file mode 100644
index 0000000000000000000000000000000000000000..1612cceda7b82dab95941f158cada99b5a5ffe96
--- /dev/null
+++ b/simc/shared/musc_with_tail.f
@@ -0,0 +1,59 @@
+ real*8 function musc_with_tail(beta,p,rad_len)
+C+_____________________________________________________________________
+!
+! MUSC_WITH_TAILS - Improved implementation of multiple scattering.
+!
+! This implementation improves on the PDG (Highland) formalism by
+! By adding a more realistic "tuned" (to HRS) heavy gaussian tail
+! by Michael Paolone.
+!
+! This formalism could further improved by replacing the Gaussian
+! with a Lynch-Dahl Gaussian as suggested by PDG, and by using the
+! true Moliere tail, but this comes at a cost of significantly
+! complicating the formula (everything because explicititly (A,Z)
+! dependent.
+!
+! Sylvester Joosten (sjoosten@anl.gov, 2019)
+!
+C-_____________________________________________________________________
+
+ implicit none
+
+! Improved formalismas in Lynch-Dahl eq. 6
+ real*8 Es, epsilon
+ parameter (Es = 13.6) !MeV
+ parameter (epsilon = 0.088)
+ real*8 rad_len
+ real*8 beta, theta_sigma, p
+ real*8 nsig_max
+ parameter(nsig_max=99.0e0) !max #/sigma for gaussian ran #s.
+ real*8 gauss1
+ real*8 norm, a
+ real*8 tail_factor
+
+! This factor adjust the total extra tail width of the wide gaussian
+! --> 1.15 to 1.2 works well for the HRS spectrometers. Adjust as
+! needed.
+ parameter(tail_factor=1.2)
+
+! Start out by calculating theta_sigma as in musc_pdg
+ theta_sigma = Es/p/beta * sqrt(rad_len) * (1+epsilon*log10(rad_len/beta**2))
+
+! Empirical fit to normalize the gaussian to have an "identical" width
+! below the defined width( good to 1-2%). Good between tail_factor
+! between 1.0 --> 20.0. Not good below 1.0
+ norm = 0.155619 / (tail_factor + 0.444642) + 1.0
+ > - exp(-tail_factor * tail_factor) / 10.0
+
+! Actual algorithm
+ theta_sigma = theta_sigma / norm
+ a = gauss1(nsig_max) * theta_sigma
+ do while (abs(a) > theta_sigma)
+ theta_sigma = tail_factor * abs(a)
+ a = gauss1(nsig_max) * theta_sigma
+ end do
+
+ musc_with_tail = gauss1(nsig_max) * theta_sigma
+
+ return
+ end
diff --git a/simc/target.f b/simc/target.f
index 2dfb428ce2b5d4e24c386b96df99b276daa86cde..5d96e7b989176fa02f8718cd19a6a01291ab0e62 100644
--- a/simc/target.f
+++ b/simc/target.f
@@ -549,30 +549,21 @@ C the perfect range, but it's easier than reproducing the generated limits here
implicit none
- real*8 Es, epsilon, nsig_max
+ real*8 musc_pdg, musc_with_tail
+
+ real*8 p, beta, teff, dangles(2)
+ real*8 dangle, r, theta_sigma, nsig_max
+ real*8 Es, epsilon
parameter (Es = 13.6) !MeV
parameter (epsilon = 0.088)
- parameter (nsig_max = 3.5)
-
- real*8 p, beta, teff, dangles(2), dangle, r
- real*8 theta_sigma
- real*8 gauss1
+ parameter(nsig_max=99.0e0)
if (p.lt.25.) write(6,*)
> 'Momentum passed to target_musc should be in MeV, but p=',p
-! Compute rms value for planar scattering angle distribution, cf. PDB
-! Note teff is thickness of material, in radiation lengths.
-
-c theta_sigma = Es/p/beta * sqrt(teff) * (1+epsilon*log10(teff))
-C Better form for beta .ne. 1, from Lynch and Dahl, NIM B58 (1991) p.6-10, Eqn. 6
- theta_sigma = Es/p/beta * sqrt(teff) * (1+epsilon*log10(teff/beta**2))
-
-! Compute scattering angles in perpendicular planes.
-! Generate two Gaussian numbers BELOW nsig_max.
- dangles(1) = theta_sigma * gauss1(nsig_max)
- dangles(2) = theta_sigma * gauss1(nsig_max)
+ dangles(1) = musc_with_tail(beta, p, teff)
+ dangles(2) = musc_with_tail(beta, p, teff)
return