cxAbyxpXsec.f

Go to the documentation of this file.

!     ratio of SxA/Sxp as a funciton of Sxp and A
!     (inelastic cross section) 
!     Need table in Cosmos/Data/DPM/sigxAbysigxp
!                       ...   /QGS/ ..
!                       ...   /EPOS/ ..
!      if xp < 15, use cxAbyxpXsecOld else use table
!
      function cxabyxpxsec(xpin, Ain) result(ratio)
      implicit none
#include "Zevhnp.h"
#include "Zmanagerp.h"
      real(8),intent(in):: xpin   ! xp inelastic cross section. (mb)
                        ! x may be any. (pi,K,p, gamma, nu )
      real(8),intent(in):: Ain    ! target mass #. May be non integer.
                            ! A = 2 ~ 210
      real(8):: ratio           !   SigxA/sigxp

      logical,save:: first=.true.
      character(40)::fileDPM="$COSMOSTOP/Data/DPM/sigxAbysigxp"
      character(40)::fileQGS="$COSMOSTOP/Data/QGS/sigxAbysigxp"
      character(40)::fileEPOS="$COSMOSTOP/Data/EPOS/sigxAbysigxp"
      character(40)::file
      integer::icon

      real(4),allocatable,save::XSratio(:,:)
      real(4),allocatable,save::XSa(:)
      real(4),allocatable,save::Aa(:)
      real(4)::xp, A, ratioS,  error
      real(4):: XSv, Av, XSr, Ar

      integer,save::iXS, iA      
      integer::iXSc, iAc
      real(8):: cxAbyxpXsecOld


      real(8):: xsxa
      if( first ) then
       !  read table sigxAbysigxp 
         if( sxabysxpopt == 1 ) then
            file =fileqgs
         elseif(sxabysxpopt == 2) then
            file =filedpm
         elseif(sxabysxpopt == 3) then
            file = fileepos
         elseif(sxabysxpopt == 4) then
           !  use old one
            goto 200
         elseif( sxabysxpopt == 0) then
            ! use default
            file =fileqgs
         else
            write(0,*) ' SxAbySxpOpt =',sxabysxpopt,' invalid'
            write(0,*) ' it must be 1 ~4'
            stop
         endif
         call copenf(tempdev, file, icon)
         if(icon /= 0) then
            write(0,*) file, " could not be opened "
            stop
         endif
         xsv =0.
         av = 0.
         ixs = 0
         ia = 0
         do while(.true.)
            read(tempdev,*, end=100) xsr, ar 
            if( xsr > xsv ) then
               ixs = ixs + 1
               xsv = xsr
            endif
            if( ar > av ) then
               ia = ia + 1 
               av = ar
            endif
         enddo
 100     continue
         rewind tempdev
         allocate( xsratio(ixs, ia) )
         allocate( xsa(ixs) )
         allocate( aa(ia) )
         do ixsc = 1, ixs
            do iac = 1, ia
               read(tempdev,*)
     *              xsa(ixsc), aa(iac), xsratio(ixsc, iac)
            enddo
         enddo
         close(tempdev)
 200     continue
         first = .false.
      endif
      if( ain == 1.0 ) then
         ratio = 1.
         return  !******************
      endif

      if( xpin < 15.d0 .or. sxabysxpopt == 4 ) then
         ratio = cxabyxpxsecold(xpin, ain)
      else
         xp = xpin
         a = ain
         call  kpolintp2s(xsa, 1, 0., aa, 1, 0.,
     *   xsratio, ixs, ixs, ia,  2, 2, xp, a, ratios, error)
         ratio = ratios
         if( ratio < 1. ) then
            ratio = cxabyxpxsecold(xpin, ain)
         endif
      endif
      end   function cxabyxpxsec

      function cxabyxpxsecold(xpin, a) result(ratio)
!       This give old value; for large xpin and small Ain
!       ration seems to be little bit smaller than other
!       code. 
      implicit none
      real(8),intent(in):: xpin ! see cxAbyxpXsec
      real(8),intent(in):: a ! see cxAbyxpXsec
      real(8):: ratio   !  xA/xp cross-section ratio
      real(8):: cinelCosByPdg

      real ca, cb, cc
      if(xpin .lt. 15.) then
         ca = .93812*a**1.0215
         cb = .65000e-02* a**0.5 + (((-.37436e-10* a +.18937e-07)*a 
     *      -.33424e-05)* a + .21646e-03)*a -.31385e-02
         ratio = ca/(1.0 + cb* xpin)
      elseif(xpin .lt. 35.) then
         cb = (0.17859e-05*a-.48592)*a +  .66909e-02  +
     *        0.48968 * a**0.99875       
         ca = (0.67395e-03*a +   .93859 )*a + .26890  +
     *        .34987e-01 * a**1.4525   
         cc = 0.925
         ratio =  ca/(1.0 + cb* xpin**cc)
      else
         ca = 80.000 * a**0.17829  + (((-.52469e-07*a +.30395e-04)*a
     *        -.69823e-02)*a +1.0517)*a +16.427
         if(a .lt. 20.) then
            cb  = 90.029* a**(-1.0219) + 
     *         (.94351e-02*a -.45131)*a+6.4772
         else
            cb = 130.79*a**(-1.0) + 
     *        (-.89803e-05*a +.52418e-02)*a -1.2073
         endif
         ratio = ca* log(1.0 + xpin/cb)/xpin
        endif
        ratio = ratio/cinelcosbypdg(a)
        ratio =  max(1.d0, ratio)
        end