cxAbyxpXsec.f File Reference

Go to the source code of this file.

Functions/Subroutines
real(8) function	cxabyxpxsec (xpin, Ain)
real(8) function	cxabyxpxsecold (xpin, a)

Function/Subroutine Documentation

cxabyxpxsec()

real(8) function cxabyxpxsec	(	real(8), intent(in)	xpin,
		real(8), intent(in)	Ain
	)

Definition at line 9 of file cxAbyxpXsec.f.

References copenf(), d0, false, kpolintp2s(), and true.

      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

Here is the call graph for this function:

cxabyxpxsecold()

real(8) function cxabyxpxsecold	(	real(8), intent(in)	xpin,
		real(8), intent(in)	a
	)

Definition at line 111 of file cxAbyxpXsec.f.

References d0, and e.

!       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)