cGetXsec.f File Reference

Go to the source code of this file.

Functions/Subroutines
subroutine	cgetxsec2 (model, pj, media, xs, mfp)
subroutine	cgetxsec (modelin, pj, media, xs, mfp)
subroutine	cxsjam (pj, media)
subroutine	cxsphits (pj, media)
subroutine	cxsdpmjet3 (pj, media)
subroutine	cxsqgsjet2 (pj, media)
subroutine	cxsepos (pj, media)
subroutine	cxssibyll (pj, media)
subroutine	cxsgheisha (pj, media)
subroutine	cxsincdpm3 (pj, media)
subroutine	cxsother (pj, media)
subroutine	cfixtarget (media)
subroutine	cworkaround (A, Z, xs, nelem)
subroutine	cfixtargetmuni (media)
subroutine	cgetcaprate (media)
subroutine	cgetphotopxs (model, pj, media, xs, mfp)

Function/Subroutine Documentation

cfixtarget()

subroutine cfixtarget

(

type(xsmedia), intent(inout)

media

)

Definition at line 498 of file cGetXsec.f.

References cerrormsg(), cworkaround(), and rndc().

Referenced by cfixtargetmuni(), cintemuon(), and csampintl().

      use modxsecmedia, xmedia=>media
      implicit none
#include "Zevhnp.h"

!
!       Fix the target element, 
!  
      type(xsmedia),intent(inout):: media  ! input/output
                        !   xsmedia information.


      real*8 u, csigma
      integer  j, ia


      if(  media%xs == smallxs .or.  media%xs == largexs ) then
         j = 1
         media%elem(j)%nsigma = media%xs
      elseif( media%noOfElem .eq. 1 ) then
         j = 1
      else
         call rndc(u)
         u = u* media%xs
         csigma = 0.
         do j = 1, media%noOfElem
            csigma = csigma + media%elem(j)%nsigma
            if(u <= csigma)  goto 10
         enddo
         write(0,*) 'media name=', media%name
         write(0,*) 'media%xs=',media%xs
         write(0,*) 'media%noOfElem=', media%noOfElem
         write(0,*) 'media%elem(:)%nsigma=',
     *          media%elem(1:j)%nsigma
         write(0,*) ' u=',u, ' csigma=',csigma, ' j=',j
         call cerrormsg('should not come here; cfixTarget',0)
      endif
 10   continue
      colelemno = j
!          int value is taken.
!c      if(model .eq. "dpmjet3" ) then
      ia = media%elem(j)%A  + 0.5
!c      else
!c         call rndc(u)
!c         ia = media%elem(j).A
!c         if(u .lt. media%elem(j).A -ia ) then
!c            ia = ia + 1
!c         endif
!c      endif
         
!      media%colZ = media%elem(j)%Z
!      media%colA = ia
!      media%colXs = media%elem(j)%nsigma/media%elem(j)%No
      targetprotonno = media%elem(j)%Z
!      TargetNucleonNo = media%elem(j)%A
      targetnucleonno = ia
      targetxs =  media%elem(j)%nsigma / media%elem(j)%No
!       Nxt needs not be called if !MacIFC  
!      (see chAcol.f and cheavyInt.f)
#if defined (MacIFC)      
      call  cworkaround(targetnucleonno, targetprotonno, targetxs,
     *   colelemno)
#endif

Here is the call graph for this function:

Here is the caller graph for this function:

cfixtargetmuni()

subroutine cfixtargetmuni

(

type(xsmedia), intent(inout)

media

)

Definition at line 585 of file cGetXsec.f.

References cfixtarget().

Referenced by csampintl().

!        fix target for muon nuclear interaction
!      In  the case of muon n.i,  x-section for each
!      element is not computed and elem(:)%nsigma is
!      not ready. So we roughly compute it and fix
!      the  target 
      use modxsecmedia, xmedia=>media
      implicit none
#include "Zevhnp.h"

      type(xsmedia),intent(inout):: media  ! input/output
                        !   xsmedia information.


      real*8 u, csigma, s0
      integer  j, ia

      s0=media%xs/sum( media%elem(:)%No * media%elem(:)%A )
      media%elem(:)%nsigma =
     *     s0 *media%elem(:)%No * media%elem(:)%A
      call cfixtarget(media)

Here is the call graph for this function:

Here is the caller graph for this function:

cgetcaprate()

subroutine cgetcaprate

(

type(xsmedia), intent(inout)

media

)

Definition at line 608 of file cGetXsec.f.

References cmucap().

Referenced by cintemuon().

      use modxsecmedia, xmedia=>media
      implicit none
#include "Zptcl.h"
#include "Zcode.h"
#include "Zevhnp.h"
          
      type(xsmedia),intent(inout):: media    ! input
 
      integer i
      real(8):: sumns, capr

      sumns = 0.
      do i = 1, media%noOfElem
         call cmucap( int(media%elem(i)%A), int(media%elem(i)%Z), 
     *            capr)
         media%elem(i)%nsigma = capr*media%elem(i)%No
         sumns = sumns + media%elem(i)%nsigma
      enddo
      media%xs = sumns   ! this is not mb x-sec. but is used
                 ! to fix the target (with media%elem(i)%nsigma

Here is the call graph for this function:

Here is the caller graph for this function:

cgetphotopxs()

subroutine cgetphotopxs	(	character*16	model,
		type(ptcl)	pj,
		type(xsmedia), intent(inout)	media,
		real(8), intent(out)	xs,
		real(8), intent(out)	mfp
	)

Definition at line 631 of file cGetXsec.f.

References cgpxsec(), and d0.

Referenced by csampgintl().

!        cgetxs for photo-hadron production
      use modxsecmedia, xmedia=>media
      implicit none
#include "Zptcl.h"
#include "Zcode.h"
#include "Zevhnp.h"
!   
      character*16 model    ! input. Active interaction model.
            ! for x-section calclation .  at present not used
      type(ptcl)::pj      ! input.  a photon
      type(xsmedia),intent(inout):: media  ! input
      real(8),intent(out):: xs   !  obtained cross-section for the
                            !  media% in mb
                    !    if xs==smallxs, no collision
                    !       xs==largexs, instant collision
      real(8),intent(out):: mfp  !  obtained  mean free path in kg/m2

      integer i
      real(8):: sumns

      sumns = 0.
      do i = 1, media%noOfElem
         call cgpxsec(media%elem(i)%A, pj%fm%p(4), xs)
         if( xs < smallxs ) then
            xs = smallxs
         endif
         if( xs == smallxs .or. xs == largexs ) then
            sumns = xs
            exit
         else
            media%elem(i)%nsigma = xs*media%elem(i)%No
            sumns = sumns + media%elem(i)%nsigma
         endif
      enddo
      media%xs = sumns

      if(media%xs /= smallxs .and. media%xs /= largexs) then
         if( media%xs <= 0. ) then
            xs = smallxs
            mfp = largexs
         else
            xs = media%xs
            mfp =1.0d0/( media%mbtoPkgrm *media%xs)
         endif
      elseif( media%xs == smallxs ) then
         xs = smallxs
         mfp = largexs
      else
         xs = largexs
         mfp = smallxs
      endif

Here is the call graph for this function:

Here is the caller graph for this function:

cgetxsec()

subroutine cgetxsec	(	character(*), intent(in)	modelin,
		type(ptcl)	pj,
		type(xsmedia), intent(inout)	media,
		real(8), intent(out)	xs,
		real(8), intent(out)	mfp
	)

Definition at line 23 of file cGetXsec.f.

References antip, cxsdpmjet3(), cxsepos(), cxsgheisha(), cxsincdpm3(), cxsjam(), cxsother(), cxsphits(), cxsqgsjet2(), cxssibyll(), cxsspecial(), d0, kneue, kneumu, and knuc.

Referenced by cgetxsec2(), and csamphadint().

! @JAXA, xs=0 leads to 0-divide error and after 10 repeats
! stop is made.  At other systems, 0 divide
! in this case has no bad effect. 
! So change is made as follows.
!    xs == smallxs is made to xs <= smallxs.
!
!         get hadronic interaction cross-section for
!    given interaction model,  projectile and medium.
!    The media information must have been set via
!    cXsecMedia.f90
!             xmedia=>media is to avoid name
!       collision of media  in modXsecMedia and
!       media argument in the subroutine def.
      use modxsecmedia, xmedia=>media
      implicit none
#include "Zcode.h"
#include "Zptcl.h"
#include "Zevhnp.h"
      character(*),intent(in):: modelin
      type(ptcl)::pj           ! input.   projectile 
      type(xsmedia),intent(inout):: media ! media 
      real(8),intent(out):: xs  ! collision xsection in mb
                    !    if xs==smallxs, no collision
                    !       xs==largexs, instant collision
      real(8),intent(out):: mfp ! mean free path in kg/m2

      character(8):: model
      
      if( pj%code == knuc ) then
         if( pj%subcode == antip .and. pj%fm%p(4) - pj%mass <= 0.) then
            xs = largexs
            media%xs = xs
            mfp = 0. 
            return  ! *******************
         endif
      elseif( pj%code == kneue .or. pj%code == kneumu ) then
         xs = smallxs
         media%xs = xs
         mfp  = largexs
         return                 ! *******************
      endif

      model = modelin      
      if( model == "special" ) then
!     model might be given in next call on return by the user;
!     (pj and target are not suited for  giving xs so that
!     some another relevant model  might have been given there for XS calc.
         call cxsspecial(pj, media, model)
      endif
       
      if( model /= "special" ) then
        select case(model)
          case( "phits" ) 
            call cxsphits(pj, media)
         case( "dpmjet3" )
            call cxsdpmjet3(pj, media)
         case( "jam" ) 
            call cxsjam(pj, media)
         case( "qgsjet2" )
            call cxsqgsjet2(pj, media)
         case( "epos" )
            call cxsepos(pj, media)
         case( "sibyll" )
            call cxssibyll(pj, media)
         case( "gheisha" ) 
            call cxsgheisha(pj, media)
         case("incdpm3")   
            call cxsincdpm3(pj, media)
         case default
            call cxsother( pj, media)
         end select
      endif
      if(media%xs < smallxs) then
         media%xs = smallxs
      endif
      if( media%xs /= smallxs .and. media%xs /= largexs) then
         mfp = 1.0d0/( media%mbtoPkgrm * media%xs)
      elseif( media%xs <= smallxs ) then
         mfp = largexs
      else
         mfp = 0.
      endif
      xs = media%xs

Here is the call graph for this function:

Here is the caller graph for this function:

cgetxsec2()

subroutine cgetxsec2	(	character(*), intent(in)	model,
		type(ptcl)	pj,
		type(xsmedia), intent(inout)	media,
		real(8), intent(out)	xs,
		real(8), intent(out)	mfp
	)

Definition at line 2 of file cGetXsec.f.

References cgetxsec().

      ! normally this is not to be used. use cGetXsec
      use modxsecmedia, xmedia=>media
      implicit none
#include "Zcode.h"
#include "Zptcl.h"
#include "Zevhnp.h"
      character(*),intent(in):: model
      type(ptcl)::pj           ! input.   projectile 
      type(xsmedia),intent(inout):: media ! media 
      real(8),intent(out):: xs  ! collision xsection in mb
                    !    if xs==smallxs, no collision
                    !       xs==largexs, instant collision
      real(8),intent(out):: mfp ! mean free path in kg/m2

      call cgetxsec(model, pj, media, xs, mfp)
      if(xs /= largexs .and. xs /= smallxs) then
         xs = xs * media%sumNo
      endif

Here is the call graph for this function:

cworkaround()

subroutine cworkaround	(	integer, intent(in)	A,
		integer, intent(in)	Z,
		real(8), intent(in)	xs,
		integer, intent(in)	nelem
	)

Definition at line 563 of file cGetXsec.f.

Referenced by cfixtarget().

!     With MacIFC,
!     TargetNucleonNo, TargetProtonNo, TargetXs
!     colElemNo
!     cannot be recognized at link time
!     from chAcol.f cheavyInt.f so we trnasfer info.
!     thru common.  (VERY STRANGE; compiler problem)
      implicit none
      integer,intent(in):: a  ! TargetNucleonNo
      integer,intent(in):: z  ! TargetProtonNo
      real(8),intent(in):: xs ! TargetXs
      integer,intent(in):: nelem ! element # of the target elemnt
                      ! in the mediumNo
#include "Zworkaround.h"

      targetnucleonno = a
      targetprotonno  = z
      targetxs   = xs
      colelemno = nelem

Here is the caller graph for this function:

cxsdpmjet3()

subroutine cxsdpmjet3	(	type(ptcl)	pj,
		type(xsmedia), intent(inout)	media
	)

Definition at line 204 of file cGetXsec.f.

References cinelx(), ctotx(), d0, knuc, and kpion.

Referenced by cgetxsec().

      use modxsecmedia, xmedia=>media
      implicit none
#include "Zptcl.h"
#include "Zcode.h"
#include "Zevhnp.h"
      type(ptcl)::pj  ! input ptcl
      type(xsmedia),intent(inout):: media    ! input/output

      integer i,  ia
      real(8):: sumns,  xs, iar

      sumns=0. 
      do i =1, media%noOfElem
         ia = media%elem(i)%A + 0.5d0
         iar = ia
         if( pj%code >= kpion .and. pj%code <= knuc ) then
            if( pj%fm%p(4) .lt.  4.1d0 ) then ! Et is better than Ek
!              call ctotx(pj, media%elem(i)%A,  xs)  !2017/sep/20 
               call ctotx(pj, iar,  xs)
!//////////////////
!               if( pj.code == 6 .and. pj.charge == -1  .and.
!     *              (pj.fm.p(4)- pj.mass) < 1d-3 ) then
!                  write(0,*) '********* ', pj.fm.p(4)-pj.mass,
!     *             xs
!                  write(0,*) ' largexs=',largexs
!               endif
!///////////////////
            else
!     call cinelx(pj, media%elem(i)%A,  media%elem(i)%Z, xs)
               call cinelx(pj, iar,  media%elem(i)%Z, xs)
            endif
         else
            call cinelx(pj, iar, media%elem(i)%Z, xs)
!            call cinelx(pj, media%elem(i)%A, media%elem(i)%Z, xs)
         endif
         if( xs <= smallxs .or. xs == largexs ) then
            sumns = xs
            exit
         else
            media%elem(i)%nsigma = xs*media%elem(i)%No
            sumns = sumns + media%elem(i)%nsigma
         endif
      enddo
      media%xs = sumns

Here is the call graph for this function:

Here is the caller graph for this function:

cxsepos()

subroutine cxsepos	(	type(ptcl)	pj,
		type(xsmedia), intent(inout)	media
	)

Definition at line 300 of file cGetXsec.f.

References cmkptc(), d0, kgnuc, and knuc.

Referenced by cgetxsec().

      use modxsecmedia, xmedia=>media
      implicit none
#include "Zptcl.h"
#include "Zcode.h"
#include "Zevhnp.h"
      type(ptcl)::pj  ! input ptcl
      type(xsmedia),intent(inout):: media    ! input/output

      integer i, ia, iz
      real(8):: sumns,  xs, tga, u, tgz

      type(ptcl)::tg

      sumns=0. 
      do i =1, media%noOfElem
         tga = media%elem(i)%A
         tgz = media%elem(i)%Z
         ia =tga + 0.5d0
         iz =tgz
!          see qgs part         
!         ia = tga
!         call rndc(u)
!         if(u .lt.  tga - ia ) then
!            ia = min(ia + 1, 209)
!         endif

         if(ia > 1 ) then
            call cmkptc(kgnuc, ia, iz, tg)
         else
            call cmkptc(knuc, ia, iz, tg)
         endif
         tg%fm%p(1:3) = 0.
         tg%fm%p(4) = tg%mass
         call ceposinioneevent(pj, tg, xs)
         if( xs <= smallxs .or. xs == largexs ) then
            sumns = xs
            exit
         else
            media%elem(i)%nsigma = xs*media%elem(i)%No
            sumns = sumns + media%elem(i)%nsigma
         endif
      enddo
      media%xs = sumns

Here is the call graph for this function:

Here is the caller graph for this function:

cxsgheisha()

subroutine cxsgheisha	(	type(ptcl)	pj,
		type(xsmedia), intent(inout)	media
	)

Definition at line 407 of file cGetXsec.f.

References cinelx(), knuc, and kpion.

Referenced by cgetxsec().

      use modxsecmedia, xmedia=>media
      implicit none
#include "Zptcl.h"
#include "Zcode.h"
#include "Zevhnp.h"
      type(ptcl)::pj  ! input ptcl
      type(xsmedia),intent(inout):: media    ! input/output

      integer i, ia
      real(8):: sumns,  xs, tga, tgz


      sumns=0. 
      do i =1, media%noOfElem
         tga = media%elem(i)%A
         tgz = media%elem(i)%Z
         if(pj%code  >= kpion .and. pj%code <= knuc ) then
            call cxsecgheisha(pj, tga,  tgz, xs)
         else
            call cinelx(pj, tga, tgz, xs)
         endif
         if( xs <= smallxs .or. xs == largexs ) then 
            sumns = xs
            exit
         else
            media%elem(i)%nsigma = xs*media%elem(i)%No
            sumns = sumns + media%elem(i)%nsigma
         endif
      enddo
      media%xs = sumns

Here is the call graph for this function:

Here is the caller graph for this function:

cxsincdpm3()

subroutine cxsincdpm3	(	type(ptcl)	pj,
		type(xsmedia), intent(inout)	media
	)

Definition at line 440 of file cGetXsec.f.

References d0.

Referenced by cgetxsec().

      use modxsecmedia, xmedia=>media
      implicit none
#include "Zptcl.h"
#include "Zcode.h"
#include "Zevhnp.h"
          
      type(ptcl)::pj      ! input.  a particle (hadronic)
      type(xsmedia),intent(inout):: media    ! input/output

      real(8):: ek, crossint
      integer:: kinc

      ek = pj%fm%p(4)- pj%mass
      if( ek > 0.2d0 ) then
!            special for inclusive treatment.  target is always air   
!         *********************************                           
         call cccode2hcode(pj, kinc)
         media%xs = crossint(kinc, ek)
      else
         media%xs = smallxs
      endif

Here is the caller graph for this function:

cxsjam()

subroutine cxsjam	(	type(ptcl)	pj,
		type(xsmedia), intent(inout)	media
	)

Definition at line 111 of file cGetXsec.f.

References cinelx(), ctotx(), d0, knuc, and kpion.

Referenced by cgetxsec().

      use modxsecmedia, xmedia=>media
      implicit none
#include "Zptcl.h"
#include "Zcode.h"
#include "Ztrackp.h"
#include "Zevhnp.h"
      type(ptcl)::pj  ! input ptcl
      type(xsmedia),intent(inout):: media    ! input/output


      integer i, ia
      real(8):: sumns,  xs, iar
      sumns=0.

      do i = 1, media%noOfElem
         ia = media%elem(i)%A + 0.5d0
         iar =  ia
         if(pj%code >= kpion .and. pj%code <= knuc ) then
            if( jamxs == 1 ) then
               call ctotx(pj, iar,  xs)
            elseif( jamxs == 0 ) then
               call cinelx(pj, iar, media%elem(i)%Z, xs)
            else
               write(0,*)
     *          ' JamXs=',jamxs, ' not usable in cxsJam'
               stop
            endif
         else
            call cinelx(pj, iar,  media%elem(i)%Z, xs)
         endif
         if( xs <  smallxs ) then
            xs = smallxs
         endif
         if( xs == smallxs .or. xs == largexs ) then
            sumns = xs
            exit
         else
            media%elem(i)%nsigma = xs*media%elem(i)%No
            sumns = sumns + media%elem(i)%nsigma
         endif
      enddo
      media%xs = sumns

Here is the call graph for this function:

Here is the caller graph for this function:

cxsother()

subroutine cxsother	(	type(ptcl)	pj,
		type(xsmedia), intent(inout)	media
	)

Definition at line 464 of file cGetXsec.f.

References cinelx(), and d0.

Referenced by cgetxsec().

      use modxsecmedia, xmedia=>media
      implicit none
#include "Zptcl.h"
#include "Zcode.h"
#include "Zevhnp.h"
          
      type(ptcl)::pj      ! input.  a particle (hadronic)
      type(xsmedia),intent(inout):: media    ! input/output

 
      integer i, ia
      real(8):: sumns, xs, iar

      sumns = 0.
      
      do i = 1, media%noOfElem
         ia =  media%elem(i)%A + 0.5d0
         iar = ia
         call cinelx(pj, iar, media%elem(i)%Z, xs)
!//////////////
!         write(0,*) 'in othert;  xs=',xs, i,  media%noOfElem
!////////////
         if( xs <= smallxs .or. xs == largexs ) then
            sumns = xs
            exit
         else
            media%elem(i)%nsigma = xs*media%elem(i)%No
            sumns = sumns + media%elem(i)%nsigma
         endif
      enddo
      media%xs = sumns

Here is the call graph for this function:

Here is the caller graph for this function:

cxsphits()

subroutine cxsphits	(	type(ptcl)	pj,
		type(xsmedia), intent(inout)	media
	)

Definition at line 156 of file cGetXsec.f.

References antip, cinelx(), kgnuc, and knuc.

Referenced by cgetxsec().

      use modxsecmedia, xmedia=>media
      implicit none
#include "Zptcl.h"
#include "Zcode.h"
#include "Zevhnp.h"
      type(ptcl)::pj  ! input ptcl
      type(xsmedia),intent(inout):: media    ! input/output


      integer i, icon
      real(8):: sumns,  xs, u, elaxs
      integer::ka, subc, ia, iz

      sumns=0. 
      ka = pj%code
      subc = pj%subcode
      do i =1, media%noOfElem
         ia = media%elem(i)%A + 0.5
         iz = media%elem(i)%Z
         if( ( ka == knuc .and. subc /= antip )
     *        .or. ka == kgnuc) then
            call cphitsxs(pj, ia, iz, elaxs,xs, icon)
!c             no need to add. xs is already total 2010.Nov.16
!cc            xs = xs + elaxs  ! phits elaxs for heavy is 0
         else
!            if( ka >= kpion .and. ka <= knuc ) then
!               if( (pj.fm.p(4)-pj.mass)  < 10.d0 )  then ! include elastic
!                  call ctotx(pj, media%elem(i)%A,  xs)
!               else
            call cinelx(pj, media%elem(i)%A, media%elem(i)%Z, xs)

!               endif
!            else
!               call cinelx(pj, media%elem(i)%A, media%elem(i)%Z, xs)
!            endif
         endif
         if( xs <= smallxs .or. xs == largexs ) then
            sumns = xs
            exit
         else
            media%elem(i)%nsigma = xs*media%elem(i)%No
            sumns = sumns + media%elem(i)%nsigma
         endif
      enddo
      media%xs = sumns

Here is the call graph for this function:

Here is the caller graph for this function:

cxsqgsjet2()

subroutine cxsqgsjet2	(	type(ptcl)	pj,
		type(xsmedia), intent(inout)	media
	)

Definition at line 252 of file cGetXsec.f.

References cinelx(), d0, kgnuc, knuc, and kpion.

Referenced by cgetxsec().

      use modxsecmedia, xmedia=>media
      implicit none
#include "Zptcl.h"
#include "Zcode.h"
#include "Zevhnp.h"
      type(ptcl)::pj  ! input ptcl
      type(xsmedia),intent(inout):: media    ! input/output


      integer i, ia
      real(8):: sumns,  xs, tga, u, tgz


      sumns=0. 
      do i =1, media%noOfElem
         tga = media%elem(i)%A
         tgz = media%elem(i)%Z
!         ia =tga 
!         call rndc(u)
!         if(u .lt.  tga - ia ) then
!            ia = min(ia + 1, 209)
!         endif
!     above  method is not so good. should  statistically be ok
!     For H2,  A  is 1.008 so some times ia becomes 2. If one makes
!     a table, 2 is always used. and dangerous.  To be accurate
!     media must be mixed A=1 andn A=2...
!     Now we treat such case roughly.  by round off (四捨五入)
         ia = tga + 0.5d0
         ia = min(ia, 209)
         if( (pj%code >= kpion .and. pj%code <= knuc) .or. 
     *        pj%code == kgnuc ) then
            call cxsecqgs(pj, ia, xs)               
         else
            call cinelx(pj, tga, tgz, xs)
         endif

         if( xs <= smallxs .or. xs == largexs ) then
            sumns = xs
            exit
         else
            media%elem(i)%nsigma = xs*media%elem(i)%No
            sumns = sumns + media%elem(i)%nsigma
         endif
      enddo
      media%xs = sumns

Here is the call graph for this function:

Here is the caller graph for this function:

cxssibyll()

subroutine cxssibyll	(	type(ptcl)	pj,
		type(xsmedia), intent(inout)	media
	)

Definition at line 346 of file cGetXsec.f.

References cmkptc(), d0, kgnuc, and knuc.

Referenced by cgetxsec().

      use modxsecmedia, xmedia=>media
      implicit none
#include "Zptcl.h"
#include "Zcode.h"
#include "Zevhnp.h"
      type(ptcl)::pj  ! input ptcl
      type(xsmedia),intent(inout):: media    ! input/output

      integer i, ia, iz
      real(8):: sumns,  xs, tga, u, tgz

      type(ptcl)::tg

      sumns=0. 
      do i =1, media%noOfElem
         tga = media%elem(i)%A
         tgz = media%elem(i)%Z
         iz =tgz
         ia = tga + 0.5d0
!     see QGS part
!                 ia =tga 
!         call rndc(u)
!         if(u .lt.  tga - ia ) then
!            ia = min(ia + 1, 209)
!         endif
         ia = min(ia, 209)
         if(ia > 1 ) then
            call cmkptc(kgnuc, ia, iz, tg)
         else
            call cmkptc(knuc, ia, iz, tg)
         endif
         tg%fm%p(1:3) = 0.
         tg%fm%p(4) = tg%mass
         if( media%name == "Air" )then
            if( i == 1) then 
               tg%subcode = 0   !sibyll can do for almost Air target only.
               call csibyllxs(pj, tg, xs)
            else
               xs =0.
            endif
         else
            call csibyllxs(pj, tg, xs)
            if( xs <= smallxs .or. xs == largexs ) then
               sumns = xs
               exit
            endif
         endif
         if( media%name ==  "Air") then
!           only 1 virtual element 'air' so weight should be 1
!     In the event generator, should not choose N target
!     but use "air" target.  i>1 xs=0
            media%elem(i)%nsigma = xs
         else
            media%elem(i)%nsigma = xs*media%elem(i)%No
         endif
         sumns = sumns + media%elem(i)%nsigma
      enddo
      media%xs = sumns

Here is the call graph for this function:

Here is the caller graph for this function: