cpdgXs.f

Go to the documentation of this file.

      module modpdgxs
      implicit none
      private 
!               if M below conflict in the user prog,
!               consider  alias like 
!              use  modpdgXs, Mpdg=>M
!               or 
!              use modpdgXs, only cSigmaT, csOfstu
      public :: csigmat, csofstu, m

      real(8),parameter:: pi= 3.14159265358979323846d0
      real(8),parameter:: hbarc=197.327d-3  ! GeV%Fermi
      real(8),parameter:: m = 2.15   ! GeV
      real(8),parameter:: eta1 = 0.462, eta2 = 0.55
      real(8),parameter:: lb = pi*(hbarc/m)**2*10.0d0   ! mb 
      real(8),parameter:: delta = 0.003056, lambda = 1.630

      contains

      function csofstu(p, Mp, Mt) result(ans)
!        s = (E+Mt)^2 - p**2 = 2*E*Mt + Mt**2 + Mp**2
      implicit none
      real(8),intent(in):: p   ! momentum  in GeV 
      real(8),intent(in):: Mp  ! proj. mass
      real(8),intent(in):: Mt  ! target mass
      real(8):: ans
      ans = (sqrt(p**2+mp**2)*2 + mt)*mt + mp**2
      end      function csofstu

      function csigmat(a,b,p) result(ans) 
!         total cross-section by pdg (ed. 2012)
      implicit none
#include "Zmass.h"
      character(*),intent(in):: a  ! one of "p", "pb", "pi+","pi-"
                                   ! "pi0"
                              !      "K+","K-" "K0" "Sig", "g"
      character(*),intent(in):: b  ! one of "p", "n", "d"
      real(8),intent(in)::p  !  momentum of  a in GeV
      real(8):: ans       ! total cross section in mb.

      real(8):: LBt, sm, s

      real(8):: Z, Y1, Y2

      lbt = lb
      select case(a)
      case('p', 'pb')
         select case(b)
         case('p')
            z = 34.71
            y1 = 12.72
            y2 = 7.35
            sm = (masp *2 + m)**2
            s = csofstu(p,masp, masp)
         case('n')
            z = 35.00
            y1 = 12.19
            y2 = 6.62
            sm = (masp + masn + m)**2
            s = csofstu(p, masp, masn)
         case('d')   
            z = 65.02
            y1= 29.04
            y2= 14.9
            sm = (masp + masd + m)**2
            s = csofstu(p, masp, masd)
            lbt =lbt*lambda
         case default
            write(0,*)
     *       'b ', b, ' for csigma(a,b) is invalid'
            stop
         end select
      case('pi+', 'pi-', 'pi0')
         select case(b)
         case('p','n')
           z = 19.02
           y1 = 9.22
           y2= 1.75
           sm = (maspic + masn + m)**2
           s = csofstu(p, maspic, masp)
         case('d') 
           z = 37.06
           y1 = 18.28
           y2 = 0.34
           sm = (maspic + masd + m)**2
           s = csofstu(p, maspic, masd)
           lbt =lbt*lambda
         case default
            write(0,*)
     *       'b ', b, ' for csigma(a,b) is invalid'
            write(0,*) ' a is ',a
            stop
         end select  
      case('K+', 'K-', 'K0')
         select case(b)
         case('p')
            z = 16.55
            y1 =  4.02
            y2 = 3.39
            sm = (maskc + masp + m)**2
            s = csofstu(p, maskc, masp)
         case('n')
            z = 16.49
            y1 = 3.44
            y2 = 1.82
            sm = (maskc + masn + m)**2
            s = csofstu(p, maskc, masn)
         case('d')   
            z = 32.34
            y1 = 7.33
            y2 = 5.59
            sm = (maskc + masd + m)**2
            s = csofstu(p, maskc, masd)
            lbt =lbt*lambda
         case default
            write(0,*)
     *       'b ', b, ' for csigma(a,b) is invalid'
            write(0,*) ' a is ',a
            stop
         end select  
      case('Sig')
         select case(b)
         case('p')
            z = 34.9
            y1 = -55.0
            y2 = -57.0
            sm = (massigmam + masp + m)**2
            s = csofstu(p, massigmam, masp)
         case default
            write(0,*)
     *       'b ', b, ' for csigma(a,b) is invalid'
            write(0,*) ' a is ',a
            stop
         end select
      case('g')
         select case(b)
         case('p','n')
            z = 34.71
            y1 = 0.0128
            y2= 0.
            sm = (masp+ masp + m)**2
            s = csofstu(p, masp, masp)
            ans = (z + lbt* (log(s/sm))**2)* delta
            sm = (masp + m)**2
            s = csofstu(p, 0.d0, masp)
            ans = ans + y1*(sm/s)**eta1
            return
         case('d')
            z = 65.02
            y1 = 0.0128
            y2= 0.
            sm = (masp + masd + m)**2
            s = csofstu(p, masp, masd)
            ans = (z + lambda* lbt* (log(s/sm))**2 )* delta
            sm = ( masd + m)**2
            s = csofstu(p, 0.d0, masd)
            ans = ans + y1*(sm/s)**eta1
            return
         case('g')
            z = 34.71
            y1 = -0.034d-4
            y2 = 0.
            sm = (masp + masp + m)**2
            s = csofstu(p, masp, masp)
            ans = (z +  lbt* (log(s/sm))**2 )* delta**2
            sm = ( masrho + m)**2
            s = csofstu(p, 0.d0, masrho)
            ans = ans + y1*(sm/s)**eta1
            return
         case default
            write(0,*)
     *       'b ', b, ' for csigma(a,b) is invalid'
            write(0,*) ' a is ',a
            stop
         end select
      case default
         write(0,*)
     *        'a ', a, ' for csigma(a,b) is invalid'
            write(0,*) ' b is ',b
            stop
      end select
      ans = z + lbt* (log(s/sm))**2 + y1*(sm/s)**eta1
      if( a == "pb" .or. a =="pi-" .or. a == "K-") then
         ans = ans + y2*(sm/s)**eta2
      else
         ans = ans - y2*(sm/s)**eta2
      endif
      end   function csigmat
      end   module modpdgxs

!      program  main
!      use modpdgXs
!      implicit none
!#include "Zmass.h"
!      character(4):: a = "g"
!      character(4):: b = "d"
!      real(8):: p,s
!      integer:: i
!      p = 5.
!      do while (p < 1.e8 ) 
!         s = csOfstu(p, 0.d0, masd)
!         write(*,*) p, cSigmaT(a, b, p), sqrt(s)
!         p = p*10.**0.1
!      enddo
!      end