csampGIntL.f

Go to the documentation of this file.

!        csampGIntL:  sample photon interaction length.
       subroutine csampgintl
       use modxsecmedia
!
!        pair, compton, photo-electric eff. photo-hadron
!
       implicit none
#include  "Zglobalc.h"
#include  "Ztrack.h"
#include  "Ztrackp.h"
#include  "Ztrackv.h"
#include  "Zevhnv.h"
#include  "Zelemagp.h"

!

       real*8  t, tkgpm2, cxai, den, cvh2den, eg
       real*8  prob, u, cmPairMFP, mfp, xs
       real*8  tcomp, tphot, tpair, tgp, tcoh
       real*8 xprob(5), txray(5)
            
       real*4  xsec(5)  !  coh, incoh,  P%E  1/(g/cm2)  
       real*8 xaimin/0.1/
!////////////////
       integer,save::nzero=0
!/////////////
!     
       eg = trackbefmove%p%fm%p(4)
       if(lpmeffect .and. eg .gt. lpmpairemin) then
!          den = cthick2den(TrackBefMove.pos.depth)
          den = cvh2den(trackbefmove%pos%height)  ! better
          call cpairlpmxsec(eg, den, prob)
          call rndc(u)
          t = -log(u) /prob
       elseif( eg .gt. 1.022d-3 ) then
          if(eg .le. excom2 ) then
             call cxrayp(eg, 4, 5,  xprob, txray)
             prob= xprob(4)+xprob(5)  !  1/r%l
             if( prob > 0. ) then
                call rndc(u)
                t = -log(u)/prob !  r%l
             else
                call cpairpath(eg, t)
             endif
          else             
             call cpairpath(eg, t)
          endif
       else
          t = infty
       endif
       call csetintinf(t*x0, .false., 'pair')

       if(eg .le. excom1 )  then
!           use Xcom data
          call cxrayp(eg,  1, 3,  xprob, txray)
          call csetintinf(txray(2)*x0, .false., 'compt')
          call csetintinf(txray(3)*x0, .false., 'photoe')
          call csetintinf(txray(1)*x0, .false., 'cohs')
       else
          if( eg .lt. maxcomptone ) then
             call ccomptpath(eg, prob, t) 
          else
             t = infty
          endif
          call csetintinf(t*x0, .false., 'compt')
          if( eg .lt. maxphotoe ) then
             call cphotoeepath(eg, prob, t)
             call csetintinf(t*x0, .false., 'photoe')
          endif
       endif
!           photo hadron prod.
       if( eg .gt. minphotoprode .and. howphotop > 0 ) then
          call cgetphotopxs(activemdl2, trackbefmove%p,
     *       media(1), xs, mfp)    ! xs in mb. mfp in kg/m2
          call rndc(u)
          tkgpm2 = -mfp*log(u)
!!!          call cghPath(eg,    tkgpm2)  ! old interface
          call csetintinf(tkgpm2, .false., 'photop')
       endif

!           magnetic  pair
       if(magpair .eq. 1 .and. eg .gt.  magpairemin) then
          xai = cxai(trackbefmove%p, mag)
          if(xai .gt. xaimin) then
             call rndc(u)
             t = -log(u) * cmpairmfp(eg, xai)
             call csetintinf(t, .true., 'mpair')
          endif
       endif
       end