cmudEdx.f

Go to the documentation of this file.

      subroutine cmudedx(flagN, flagBr, flagPr,  Emu, dEdx, dEdxF)
      implicit none

!          gives muon sum of the energy loss due to direct pair creation,
!          bremsstrahlung and nuclear interaction.
!          dE/dx  by ionization is not included here.
!
      integer flagN  ! input.  specifies how to treate nuclear interaction.
                     !         0--> n.i is neglected always.
                     !         1-->int(v0:vmax) of Emu* v*dsigma/dv is
                     !             put in dEdx.v0~0, vmax~1.0.
                     !             (v=E_loss/Emu)
                     !            I.e., n.i is treated as a continuous process.
                     !         2/3-->int(0:vmin) of Emu* v*dsigma/dv is
                     !             put in dEdx. vmin=cnst.muNVmin(=10^-3)
                     !         (discrete process  by  v>vmin is sampled
                     !          for each case; 2--> no n.i is followed.
                     !          3--> n.i is followed)
                     !          
      integer flagBr  ! input.  The same meaning for bremsstrahlung as flagN.
      integer flagPr  ! input   //                   direct pair creation.
      real*8  Emu     ! input. muon total energy
      real*8  dEdx    ! output.   dE/dx GeV/(g/cm^2). sum of loss due to
                      !           the three  process.  restricted one 
                      !           if frag is 2/3
      real*8  dEdxF !  output  Full dE/dx. if flag =1, same as dEdx
!
      real*8  dEdx1, dEdx2,  dEdx3
      real*8  dEdx1F, dEdx2F,  dEdx3F
!
      call cmundedx(flagn,  emu, dedx1)
      if(flagn .ne. 1) then
         call cmundedx(1,  emu, dedx1f)
      else
         dedx1f = dedx1
      endif

      call cmubrdedx(flagbr, emu, dedx2)
      if(flagbr .ne. 1) then
         call cmubrdedx(1, emu, dedx2f)
      else
         dedx2f = dedx2
      endif
      call cmuprdedx(flagpr, emu, dedx3)
      if(flagpr .ne. 1) then
         call cmuprdedx(1, emu, dedx3f)
      else
         dedx3f = dedx3
      endif
      dedx = dedx1 + dedx2 + dedx3
      dedxf = dedx1f + dedx2f + dedx3f
      end
!     ********************
      subroutine cmundedx(flag, Emu, dEdx)
!     ********************
      implicit none
#include "Zcmuint.h"
!         compute dE/dx of muon by nuclear interaction
      integer flag         ! input. 0--> dEdx is always 0.
                           !        1--> dEdx is for all v
                           !        2--> dEdx is for v<vmin
                           !             nuclear interaction is 
                           !             not followed
                           !        3--> dEdx is for v<vmin
                           !             n.i is followed
      real*8  Emu          ! input. muon total energy in GeV.
      real*8  dEdx         ! output. dE/dx GeV/(g/cm2)


      real*8 ale, pw

      if(flag .eq.0 .or.  (emu .lt. munemin)) then
         dedx = 0.
      else
         if(emu .gt.  munemax1) then
            ale = log10(munemax1)
         else
            ale =  log10(emu)
         endif
         if(flag .eq.  1) then
!              all is  continuous loss
            call kintp3(mundedxt, 1, muntxt, munlemin,
     *       mundetx, ale, dedx)
            pw = munpwdedxt
         elseif(flag .eq. 2 .or. flag .eq. 3) then
!              v<vmin is regarded as continuous loss
            call kintp3(mundedx0, 1, muntxt, munlemin,
     *      mundetx, ale, dedx)
            pw = munpwdedx0
         else
            call cerrormsg('flag is invalid for cmuNdEdx',0)
         endif
         if(emu .gt. munemax1) then
            dedx = dedx*(emu/munemax1)**pw
         endif
         dedx = dedx * emu
      endif
      end
!     *******************
      subroutine cmubrdedx(flag, Emu, dEdx)
!     *******************
      implicit none
#include "Zcmuint.h"
!         compute dE/dx of muon by brems
      integer flag         ! input. 0--> dEdx is always 0.
                           !        1--> dEdx is for all v
                           !        2/3--> dEdx is for v<vmin
      real*8  Emu          ! input. muon total energy in GeV.
      real*8  dEdx         ! output. dE/dx GeV/(g/cm2)


      real*8 ale

      if(flag .eq.0 .or.  (emu .lt. mubremin)) then
         dedx = 0.
      else
         if(emu .gt.  mubremax1) then
            ale = log10(mubremax1)
         else
            ale =  log10(emu)
         endif
         if(flag .eq.  1) then
!              all is  continuous loss
            call kintp3(mubrdedxt, 1, mubrtxt, mubrlemin,
     *      mubrdetx, ale, dedx)
         elseif(flag .eq. 2 .or. flag .eq. 3) then
!              v<vmin is regarded as continuous loss
            call kintp3(mubrdedx0, 1, mubrtxt, mubrlemin,
     *      mubrdetx, ale, dedx)
         else
            call cerrormsg('flag is invalid for cmuBrdEdx',0)
         endif
         dedx = dedx * emu
      endif
      end

!     *******************
      subroutine cmuprdedx(flag, Emu, dEdx)
!     *******************
      implicit none
#include "Zcmuint.h"
!         compute dE/dx of muon by pair creation
      integer flag         ! input. 0--> dEdx is always 0.
                           !        1--> dEdx is for all v
                           !        2/3--> dEdx is for v<vmin
      real*8  Emu          ! input. muon total energy in GeV.
      real*8  dEdx         ! output. dE/dx GeV/(g/cm2)


      real*8 ale

      if(flag .eq.0 .or.  (emu .lt. mupremin)) then
         dedx = 0.
      else
         if(emu .gt.  mupremax1) then
            ale = log10(mupremax1)
         else
            ale = log10(emu)
         endif
         if(flag .eq.  1) then
!              all is  continuous loss
            call kintp3(muprdedxt, 1, muprtxt, muprlemin,
     *      muprdetx, ale, dedx)
         elseif(flag .eq. 2 .or. flag .eq. 3) then
!              v<vmin is regarded as continuous loss
            call kintp3(muprdedx0, 1, muprtxt, muprlemin,
     *      muprdetx, ale, dedx)
         else
            write(0,*) ' flag = ', flag
            call cerrormsg('flag is invalid for cmuPrdEdx',0)
         endif
         dedx = dedx * emu
      endif
      end