implemented better multiple scattering, will need testing

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

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

simc/shared/musc_pdg.f

+	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 

simc/shared/musc_with_tail.f

+	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 

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