cmulScat.f

Go to the documentation of this file.

!     ****************************************
!     *                                                          *
!     * cmulScat: multiple Coulomb scattering
!     *                                                          *
!     ****************************************
!
!
      subroutine cmulscat(theta)
      use modxsecmedia
      use modmcscontrol
      implicit none
!       Using  cTrack and Move, compute scattering angle 
!
#include "Zglobalc.h"
#include  "Zcode.h"
#include  "Ztrack.h"
#include  "Ztrackv.h"
#include  "Zelemagp.h"

      real*8 theta         ! output. sampled angle in radian.
!
!            
!
      integer cond
      real(8):: sgr
      integer ic

!             sample theta
      if(intinfarray(processno)%length .gt. 0.) then
         if(activemcs == "El_con") then
            !  total length in kg/m2; sample theta                   
            sgr =intinfarray(processno)%thickness
            call celsepacondense(sgr, theta)
         elseif(abs(moliere) >= 2)  then   ! abs is historical reason
            call csampmol(media(1), theta, cond)  ! rigorous Moliere
            if(cond /= 0 ) then
               call cmulscat2(theta)   ! Gauss
            endif
         elseif(abs(moliere) == 1)  then
            call cmulscat1(theta)      ! simplified Moliere
         elseif( moliere == 0 ) then
            call cmulscat2(theta)      ! Gauss
         endif
      else
         theta = 0.
      endif
      end
!     **************************
      subroutine cmulscat2(theta)
      implicit none
!       Gaussian scattring
!
#include  "Zglobalc.h"
#include  "Ztrack.h"
#include  "Ztrackv.h"
#include  "Zelemagp.h"

      

      real*8 theta ! output. sampled spatial angle in radain.

      real*8 tetarms, g1, g2, u, beta2, dt
      integer nc
      
      real*8 hpi 
      parameter(hpi = pi/2.)

      g1 = trackbefmove%p%fm%p(4)/movedtrack%p%mass
      g2 = movedtrack%p%fm%p(4)/movedtrack%p%mass
      beta2 = 1.d0 - 1.d0/g1/g2
      if(beta2 .le. 0.) then
         tetarms = 0.
      else
         dt = intinfarray(processno)%thickness/ x0  ! r%l
         if(dt .gt. 1.d-3) then
            tetarms = es/trackbefmove%p%fm%p(4)*
     *       abs(movedtrack%p%charge) *
     *       sqrt(dt/beta2)*(1.0 + 0.038*log(dt))
         else
            tetarms = es/trackbefmove%p%fm%p(4)*
     *      abs(movedtrack%p%charge) * sqrt(dt/beta2)
         endif
      endif
      theta = pi
      nc = 0
      do while(theta .gt. hpi)
         if(nc .gt. 10) then
!              tetarms seems too large
            theta = u**0.1 * hpi  ! give some value 
         else
            call rndc(u)
            theta = sqrt(-log(u))* tetarms
            nc = nc +1
         endif
      enddo
      end
      subroutine cmulscat1(theta)
      implicit none
#include  "Zcode.h"
#include  "Ztrack.h"
#include  "Ztrackv.h"
#include  "Zelemagp.h"
!           by Moliere thoery 

      real*8 theta
      
      real*8 g1, g2, leng
      integer cond

      real*8 t, tmp, avx, avy, disp, cs, sn, e1, e2, d1, d2
      real*8 rho, cvh2den
      integer chg
 
      e1 = trackbefmove%p%fm%p(4)
      g1 = e1/trackbefmove%p%mass
      e2 = movedtrack%p%fm%p(4)
      g2 = e2/movedtrack%p%mass
      rho = cvh2den(
     *      (trackbefmove%pos%height+movedtrack%pos%height)/2.d0
     *      )

      leng = intinfarray(processno)%length  ! in m

      chg =  movedtrack%p%charge
      call cmoliere(rho,  chg, movedtrack%p%mass, g1, g2,
     *     leng, theta, cond)
      if(cond .ne. 0) then
!           Moliere theory cannot be applied
         call cmulscat2(theta)
      endif
      end

      subroutine celsepacondense(sgr, theta)
      use modmcscontrol
      use modsoftmcs
      use modcmcs
      use modxsecmedia
      implicit none

      real(8),intent(in):: sgr  ! path length in kg/m2
      real(8),intent(out):: theta ! sampled angle
      
      real(8)::cm2topgrm,s0,s1, s2, ls1, ls2 
      real(8):: mu, cosa, sina
!!!      call ciniFsoft(KEeV, sgr)   is for integration
!!!    from 0 to muc.  We need here from 0 to 1. so
!!!   do below
!                         1 is to get s in cm2
      call ctpxs(tpxsnow, 1,  keev, s0, s1, s2)
!                              not 0.1 (/(kg/m2))-->(/g/cm2)
      cm2topgrm = media(1)%mbtoPkgrm/ 1d-27 
      ls1 =1.0/(s1*cm2topgrm)
      ls2 =1.0/(s2*cm2topgrm)

      avemu = (1.0d0 - exp(-sgr/ls1))/2
      avemu2 = avemu - (1- exp(-sgr/ls2) )/6.d0
      if( avemu > 0.49999999d0 ) then
          ! mu is uniform (asoft =1 bsoft =1 is also ok)
         bsoft = 0.5d0
         asoft = 0.5d0
      else
         bsoft =( 2*avemu - 3*avemu2 )/(1-2*avemu)
         asoft = (1-2*avemu) + bsoft
      endif
      call csampsoftmcsang(mu, cosa)
      theta = acos(cosa) 
      end