csampNEPIntL.f

Go to the documentation of this file.

!     ******************************************************************
!     *                                                                *
!     * csampNEPIntL: samples integraction length for a given
!     *              non e, gamma  particle in a given material air
!     *                                                                *
!     ******************************************************************
!c            From Epicsv7.0. this is not used from Epcis
!       I.e. this is for Air
!     *****************************************************************
      subroutine  csampnepintl
      use modxsecmedia
      implicit none
#include  "Zglobalc.h"
#include  "Zcode.h"
#include  "Ztrack.h"
#include  "Ztrackv.h"
#include  "Zevhnp.h"

!     **************************************************
!
      real*8 mfp, xs, length

!            
!        m.f.p (kg/m**2) = abn /xsec(mb)
#if defined (LOOPCHK)
      integer ksave/0/, esave/-1./, ncount/0/
      integer ka
      save
      ka = trackbefmove%p%code

      if(ksave .eq. ka .and. esave .eq.  
     *   trackbefmove%p%fm%p(4)) then
         ncount = ncount +1
         if(ncount .gt. 10) then
            write(0,*)'  ncount=',ncount, 
     *        ' ka=',ka, ' e=',esave, ' mass=',
     *      trackbefmove%p%mass
         endif
      else
         ksave =ka
         esave = trackbefmove%p%fm%p(4)
         ncount =0
      endif
!///////////////////
#endif

      call cdecayleng(trackbefmove, length)

      if(length .ne. infty) then
         call csetintinf(length, .true., 'decay')
      endif
!      nnb ddb rho phi omega ---> length=0 (instant decay)
!              or stopping mu pi etc -> length=0 
      if(length .gt. 0.) then 
         if(trackbefmove%p%charge .ne. 0) then
!                  heavy (non e+/e-) knockon
!             if Freec and mu or had makes first hadronic
!             interaction, this call will be non-effective
            call cknockonh
         endif
!                  non stopping-decay paticle
         if( trackbefmove%p%code .eq. kmuon) then
!                   muon; pair, brems, n.i   
            call csampmuint(xs, mfp)
         else   
!                hadronic  collisions
            call csamphadint(xs, mfp)

!           fixing target will be done after collision
!          is really selected.
!            if(xs >  smallxs .and. xs .ne. largexs) then
!                 fix target
!               call cfixTarget(media(1))
!            elseif( xs == largexs ) then
!               call cfixTarget(media(1))   !  use big xs 
                         ! to fix the target. pbar annihilation
                         ! probably A independent.  
!            endif
         endif
      endif
      end

      subroutine csampmuint(xs, mfp)
      use modxsecmedia
      implicit none
#include  "Zglobalc.h"
! #include  "Zcode.h"
#include  "ZbpCnst.h"
#include  "Ztrackp.h"
#include  "Ztrack.h"
#include  "Ztrackv.h"
! #include  "Zheavyp.h"
#include  "Zelemagp.h"
#include  "Zevhnp.h"
! #include  "Zevhnv.h"
#include  "Zcmuint.h"
      real*8   xs, mfp, u, csigma, s0

      real*8  prob, path, collkgram

      integer::i

      xs = smallxs
      if(mupr .ge. 2 .and.
     *     trackbefmove%p%fm%p(4) .gt. mupremin) then
         call cmuprsmpp(trackbefmove%p%fm%p(4), prob, path)
!          prob is /r.l  path is in r.l     
!          prob   /r.l /X0    --> prob/ (kg/m2) *10 -> prob/(g/cm2)
!          = mb*4.138e-5
!            mb2pg = 1.e-27 * Avogn/A = 4.138e-5
!          xs =prob*10/X0/mb2pg = prob*10/X0*A/Avogn*1.e27
!             = prob/X0*A/Avogn*1.e28
!             = prob/X0*A*A2deninv    
!       (   A2deninv = 1.d28/Avogn   )
!         xs =max( prob/X0* A2deninv*TargetMassN, smallxs)
!           n mfp xs =1    
         xs =max( prob/mbtopx0, smallxs)
         mfp = 1./(media(1)%mbtoPkgrm* xs)
         collkgram = path * x0  ! kg/m^2
         call csetintinf(collkgram, .false., 'mupair')
      endif
      if(mubr .ge. 2 .and. 
     *     trackbefmove%p%fm%p(4) .gt. mubremin) then
               
         call cmubrsmpp(trackbefmove%p%fm%p(4), prob, path)
         xs =max( prob/mbtopx0, smallxs)
         mfp = 1./(media(1)%mbtoPkgrm* xs)
         collkgram = path * x0 ! kg/m^2
         call csetintinf(collkgram, .false., 'mubrem')
      endif
      if(muni .ge. 2 .and. 
     *     trackbefmove%p%fm%p(4) .gt. munemin) then
         call cmunsmpp(trackbefmove%p%fm%p(4), prob, path)
         xs =max( prob/mbtopx0, smallxs)
         media(1)%xs = xs  ! for MuNI=3
         mfp = 1./(media(1)%mbtoPkgrm* xs)
         collkgram = path * x0 ! kg/m^2
         if( xs .ne. smallxs .and.
     *        .not. freec .and. zfirst%pos%depth .eq. 0. ) then
            collkgram = 0.
         endif
         call csetintinf(collkgram, .false., 'munuci')

!          next is postponed after munuci is really selected
!         if( xs > smallxs .and. xs .ne. largexs ) then
!            call cfixTargetMuNI(media(1))
!         endif
      endif
      end


      subroutine csamphadint(xs, mfp)
      use modxsecmedia
      implicit none
#include  "Zglobalc.h"
#include  "Ztrackp.h"
#include  "Ztrack.h"
#include  "Ztrackv.h"
#include  "Zevhnp.h"
#include  "Zevhnv.h"
#include  "Zcmuint.h"
      real(8),intent(out):: xs !  mb
      real(8),intent(out):: mfp !  kg/m2

      real*8  prob, path, collkgram, uta, u
      call cgetxsec(activemdl2,  trackbefmove%p, media(1),
     *     xs, mfp)
!     for hadronic interaction 
!                 we can fix the fist col. point at the input
!                 point when Freec=f
      if(xs .eq. smallxs) then
         collkgram = infty
      elseif(.not. freec .and. zfirst%pos%depth .eq. 0. ) then
!         forced collision at the inut deth.
         collkgram = 0.
      elseif( xs .eq. largexs) then
!                  say stopping pbar
         collkgram = 0.
      else
!                 sample interaction length
         call rndc(u)
         collkgram=-mfp*log(u)
      endif
      call csetintinf(collkgram, .false., 'coll')

      end


      subroutine cknockonh
      implicit none
#include  "Zglobalc.h"
#include  "Zcode.h"
#include  "Ztrack.h"
#include  "Ztrackv.h"
#include  "Zelemagp.h"
#include  "Zevhnv.h"
      real*8 prob, path
!              knock on by non e+/e- charged ptcl
      call cknockp(trackbefmove%p, prob, path) ! path in r%l
      if(prob .gt. 0.d0) then
         call csetintinf(path *x0, .false., 'knock')
      else
         call csetintinf(infty,   .false., 'knock')
      endif
      end