cdedxatHE.f File Reference

Go to the source code of this file.

Functions/Subroutines
subroutine	cdedxathe (aPtcl, rhoin, dedt)

Function/Subroutine Documentation

cdedxathe()

subroutine cdedxathe	(	type(ptcl)	aPtcl,
		real*8	rhoin,
		real*8	dedt
	)

Definition at line 23 of file cdedxatHE.f.

References betasq, d, d0, d3, parameter(), w0, and wlg0.

         implicit none
#include  "Zptcl.h"
#include  "ZdedxAir.h"
         type(ptcl):: aptcl
         real*8  rhoin, dedt
!
!

         real*8 emass, emass2
         parameter(emass=.511d0, emass2=emass**2)
!
         real*8  ek, amass, a, b,  x0, c, x1, sa, rho, rl, e, gi
         real*8  pbmc, cb, dltx, p2, wm, wlg, x
!
!          see      sternheimer's consts. by p.r.b vol.3 (1971)3681
!                 to mev unit
       if(jdef .eq. 0) then
          write(*, *) ' cdedxEleci must be called'
          stop 9999
       endif

         if(aptcl%charge .ne. 0) then
              ek=aptcl%fm%p(4) *1.d3
              amass=aptcl%mass * 1.d3
              a=7.68d-2
              b=17.87d0
              if(ek .lt. 50.d0*amass) then
                  x0=1.884d0
                  c=-10.97d0
                  x1=4.d0
                  sa=.247d0
              else
                  rho=min( max(rhoin, 1.d-7), 5.d-3)
                  rl=log10(rho)
                  x0= (((-.0165d0*rl-0.305d0)*rl-1.94d0)*rl
     *             -5.41d0)*rl-3.87d0
                  c=-4.0635d0+2.303d0*rl
                  if(rho .gt. 2.78605d-4) then
                      x1=4.d0
                      sa=((-.6872064d-01*rl-.5340530d0 )*rl
     *               -1.521159d0 )*rl-1.365158d0
                  else
                      x1=5.d0
                      sa=(((.04256d0*rl+0.7888d0)*rl+5.465d0)*rl
     *                   + 16.642d0)*rl+    18.855d0
                  endif
              endif
!                     energy in mev unit
              e=(ek+amass)
              gi=amass/e
              betasq = 1.d0 - gi**2
!                      x=log10(p/mc)
              pbmc =      (e/amass)**2 - 1.d0
              if(pbmc .gt. 0.) then
                  x=log10(pbmc)/2.0d0
                  cb=-c
                  if(x .lt. x0) then
!                        4.606x - dlt
                      dltx=4.606d0*x
                  elseif(x .lt. x1) then
                      dltx=cb - (x1-x)**3 * sa
                  else
                      dltx=cb
                  endif
                  p2=e**2 -amass**2
                  wm=2*emass*p2/( amass**2+ emass2+ emass*e*2)
!
                  if(wm .gt. w0) then
                     wlg = wlg0
                  else
                     wlg=log(wm)
                  endif
                  dedt=a/betasq *( b+0.693d0+wlg -2.d0*betasq + dltx)
!                                convert it to GeV/(kg/m2)
                  dedt=dedt *1.d-4 *aptcl%charge**2
            else
                  dedt=1.d-9
            endif
         else
            dedt = 0.
         endif

Here is the call graph for this function: