cheavyInt.f File Reference

Go to the source code of this file.

Functions/Subroutines
subroutine	cheavyint (pj, ia, iz, a, ntp)
subroutine	csampfragments (pj, ia,
subroutine	cresetnucchg (nuc, noOfNuc, z)
subroutine	csampintenuc (pj, ia, noOfNuc, noOfInteNuc)
subroutine	csampfragmom (a, nf)
subroutine	csampfragpt (aPtcl, pt)
subroutine	cmakefragtbl
subroutine	chg2massn (hg, massnum)
subroutine	chg2charge (hg, charge)
subroutine	ccode2hvgrp (code, hg)
subroutine	ccode2mass (code, mass)

Function/Subroutine Documentation

ccode2hvgrp()

subroutine ccode2hvgrp	(	integer	code,
		integer	hg
	)

Definition at line 549 of file cheavyInt.f.

References cerrormsg(), kalfa, and kiron.

      implicit none
#include "Zcode.h"
#include "Zheavyp.h"
      integer code, hg
      if(code .ge. kalfa .and. code .le. kiron) then
         hg = code2heavyg(code)
      else
         call cerrormsg(
     *   'ccode2hvgrp should not be used for code # He ~Fe',0)
      endif

Here is the call graph for this function:

ccode2mass()

subroutine ccode2mass	(	integer	code,
		real*8	mass
	)

Definition at line 561 of file cheavyInt.f.

References cerrormsg(), kalfa, kiron, masn, and masp.

      implicit none
#include "Zcode.h"
#include "Zheavyp.h"
#include "Zmass.h"
      integer code, hg
      real*8 mass
      if(code .ge. kalfa .and. code .le. kiron) then
         hg = code2heavyg(code)
         mass = heavyg2massn(hg) * (masp+masn)/2
      else
         call cerrormsg(
     *   'ccode2mass should not be used for code # He ~Fe',0)
      endif

Here is the call graph for this function:

cheavyint()

subroutine cheavyint	(	type(ptcl)	pj,
		integer	ia,
		integer	iz,
		type(ptcl), dimension(*)	a,
		integer	ntp
	)

Definition at line 16 of file cheavyInt.f.

References kgnuc, and maxheavymassn.

Referenced by cinteheavy().

#if !defined (MacIFC)
        use modxsecmedia
#endif
        implicit none

#include  "Zcode.h"
#include  "Zptcl.h"
#include  "Zcoord.h"
#include  "Zheavyv.h"
#include  "Zevhnp.h"
#include  "Zevhnv.h"
#include  "Zair.h"

!///////////
!        logical deb
!        common/cdebug/ deb
!//////////
!
        external cblkheavy
!
!
!          ia: mass no. of the target
!          iz: charge of the target

        integer ia, iz, ntp
        type(ptcl):: pj, a(*)
!
!        type(coord):: dir
!        real*8  dummylen
!
!        to store heavy fragment            nucleon
        type(ptcl):: frga(maxheavymassn), nuc(maxheavymassn)

        type(ptcl):: fraga(maxheavymassn), intnuca(maxheavymassn),
     *    nonintnuca(maxheavymassn)
!

        integer nooffragments, noofnuc, noofintenuc
        save  nooffragments, noofnuc, noofintenuc

        integer nooffrag, noofintn, noofnonintn

        integer  i, n, j
!///////////////
!        integer cumev, cev
!////////////
!         /////
!         real*8 sumerg

!         call ctestOnShell(' heavy bef frag', pj)
!         sumerg = 0.
!         ////

!     
! MacIFC  ld error  ; also for TargetXs here and chAcol.f)
!  (No error when complile ephook.f but for
!   testCnf1.mk error occurs)
!////        integer:: dummy
!///         dummy = TargetNucleonNo
!Undefined symbols for architecture x86_64:
!  "_modxsecmedia_mp_targetnucleonno_", referenced from:
!      _cheavyint_ in libcosmos.a(cheavyInt.o)
!  "_modxsecmedia_mp_targetxs_", referenced from:
!      _cheavyint_ in libcosmos.a(cheavyInt.o)
!      _chacolx_ in libcosmos.a(chAcol.o)
!ld: symbol(s) not found for architecture x86_64
!c//////////////////
#if defined (MacIFC)
#include "Zworkaround.h"
#endif

        ntp=0
!
!            ** fragmentation with some nucleon interaction **
!                 sample and set fragmentation ptcls
!      noOfFragments:  heavy fragments
!      noOfNuc:  all nucleons
!      noOfInteNuc: interacting nucleons.
!
!                                          
!        if(IntModel .eq. 'int2' .and.
!     *      ( pj.code .eq. ktriton .or. pj.code .eq. kdeut) ) then

        if( activemdl .eq. 'dpmjet3') then
           if((pj%fm%p(4) - pj%mass)/pj%subcode .gt. 5.1) then
              call cdpmjet(pj, ia, iz, a, ntp)
           else
!               rescue
              call cjamevent(pj, ia, iz, targetxs, a, ntp)
           endif
        elseif( activemdl .eq. 'jam') then
           call cjamevent(pj, ia, iz, targetxs, a, ntp)
        elseif( activemdl .eq. 'phits') then
           call cphits(pj, ia, iz, targetxs, a, ntp)
        elseif(activemdl .eq. 'qgsjet2') then
           call cqgsjet(pj, ia, iz, a, ntp)
        elseif(activemdl .eq. 'epos') then
           call ceposgenoneevent(pj, ia, iz, a, ntp)
        elseif(activemdl .eq. 'sibyll') then
           call csibyllevent(pj, ia, iz,a, ntp)
        elseif(activemdl .eq. 'gheisha' .and.
     *        pj%code .eq. kgnuc  .and. (pj%subcode .eq. 2 .or.
     *        pj%subcode .eq. 3) ) then
!                  Gheisha should be called directly
           call chagheisha(pj, ia, iz, a, ntp)
           nooffragments = 0
           noofnuc = 0
           noofintenuc = 0
        else
!              resque
           call cjamevent(pj, ia, iz, targetxs, a, ntp)           
!!           call csampFragments(pj, ia, frga, noOfFragments,
!!     *       nuc, noOfNuc, noOfInteNuc)
!!           if(SkipPtclGen .eq. 0) then
!!c                 interaction of interating nuc.
!!              do  i = 1, noOfInteNuc
!!                 j=i + noOfNuc - noOfInteNuc
!!                 call chAcol(nuc(j), ia, iz, a(ntp+1), n)
!!                 ntp=ntp+n
!!              enddo
!!c               all business so far is done in the frame where
!!c               z axis is the incident
!!c               make rotation
!!              call crot3mom( pj, a, ntp )
!!           endif
!!           call crot3mom( pj, nuc, noOfNuc) ! for all nucleons. see entry below
!!
!!           call crot3mom( pj, frga, noOfFragments)
!!
!!
!!c               move non-interacting nucleon and fragments
!!           do  i = 1, noOfNuc - noOfInteNuc 
!!              a(ntp+1) = nuc(i)
!!              ntp = ntp +1
!!           enddo
!!           do  i = 1, noOfFragments
!!              a(ntp + 1) = frga(i)
!!              ntp = ntp + 1
!!           enddo
        endif
        return
!     ****************** inquire the fragments at heavy interation
      entry cqhvyintf(fraga, nooffrag)
!     ******************* 
!        note if MovedTrack.p.code is not heavy, this gives wrong result

      nooffrag = nooffragments
!         move fragments
      do i = 1, nooffrag
         fraga(i) = frga(i)
      enddo
      return
!     ******************** inquire the interacting nuc.
      entry cqhvyintin(intnuca, noofintn)
!     *******************
      noofintn = noofintenuc
!        move  interacting nucleons
      do i = 1, noofintn
         intnuca(i) = nuc(i+noofnuc - noofintenuc)
      enddo
      return
!     ******************inquire non interacting  nuc.
      entry cqhvyintnin(nonintnuca,  noofnonintn)
!     ******************
      noofnonintn = noofnuc - noofintenuc
!        move non-interacting nucleons
      do i = 1, noofnonintn
         nonintnuca(i) = nuc(i)
      enddo

Here is the caller graph for this function:

chg2charge()

subroutine chg2charge	(	integer	hg,
		integer	charge
	)

Definition at line 542 of file cheavyInt.f.

      implicit none
#include "Zcode.h"
#include "Zheavyp.h"
      integer charge, hg
      charge = heavyg2charge(hg)

chg2massn()

subroutine chg2massn	(	integer	hg,
		integer	massnum
	)

Definition at line 535 of file cheavyInt.f.

      implicit none
#include "Zcode.h"
#include "Zheavyp.h"
      integer massnum, hg
      massnum = heavyg2massn(hg)

cmakefragtbl()

subroutine cmakefragtbl

(

)

Definition at line 507 of file cheavyInt.f.

References maxheavyg.

Referenced by csampfragments().

      implicit none

#include  "Zcode.h"
#include  "Zheavyp.h"
#include  "Zheavyc.h"
!
      integer i, j
!
!          do below for nucleus group 1 to maxHeavyG
!
       do   i=1, maxheavyg
!              FragmentaionTbl(i,j) containes 
!              <no. of nucleus of group j> when
!              nucleus of group i fragments  (j<=i)
!                make cumulative table
          cfragmentationtbl(i, 1) = fragmenttbl(i, 1)
          do   j=1,i-1
             cfragmentationtbl(i,j+1)=fragmenttbl(i,j+1)
     *          +cfragmentationtbl(i,j)
          enddo
!                 normalzie
          do  j=1, i  
             cfragmentationtbl(i, j) = cfragmentationtbl(i,j)/
     *          cfragmentationtbl(i,i)
          enddo
       enddo

Here is the caller graph for this function:

cresetnucchg()

subroutine cresetnucchg	(	type(ptcl), dimension(noofnuc)	nuc,
		integer	noOfNuc,
		integer	z
	)

Definition at line 361 of file cheavyInt.f.

Referenced by csampfragments().

!
       implicit none
!----       include '../../Zptcl.h'
#include  "Zptcl.h"
       integer z, noofnuc
       type(ptcl):: nuc(noofnuc)
!
       integer i
!
       do i = 1, z
          nuc(i)%charge = 1
       enddo
       do i = z+1, noofnuc
          nuc(i)%charge = 0
       enddo

Here is the caller graph for this function:

csampfragments()

subroutine csampfragments	(	type(ptcl)	pj,
		integer	ia
	)

Definition at line 219 of file cheavyInt.f.

References cmakefragtbl(), cmkptc(), cresetnucchg(), csampfragmom(), csampintenuc(), d0, false, kfrge(), knuc, maxheavyg, regptcl, and rndc().

     *       fra,  nooffragments, nuc, noofnuc, noofintenuc)
        implicit none

#include  "Zptcl.h"
#include  "Zcode.h"
#include  "Zheavyp.h"
#include  "Zheavyc.h"


        integer nooffragments, noofnuc, noofintenuc, ia
!         ia: # of target nucleons
        type(ptcl):: pj, fra(*), nuc(*)
        integer ihg, mno, mx, mn, msumf, ihgf, jcon,
     *          hchg
        real*8  epn, u
        logical first/.true./
        integer i, zfrag
!
      if(first) then
         call cmakefragtbl
         first = .false.
      endif

      if( (pj%fm%p(4)-pj%mass)/pj%subcode .lt. 0.1) then
         fra(1) = pj
         nooffragments = 1
         noofnuc = 0
         noofintenuc = 0
         return   !   ********
      endif
!
!         get heavy group index from the incident charge  
      ihg=charge2heavyg(pj%charge)

!          # of nucleons
!      mno=HeavyG2massN(ihg)
      mno = pj%subcode
!       z 
!      hchg=HeavyG2charge(ihg)

       hchg=pj%charge

!       t.e  energy per nucleon
      epn=pj%fm%p(4)/mno
!         margin for mass no. conservation (for first trial)
      mx = mno
      mn = mno - 3

!
!             retry if final sum of mass exceeds incident one
!      *** until loop*** 
      do while (.true.)
!            sum of mass of fragments
         msumf = 0
!            no. of heavy fragments
         nooffragments = 0
!            no. of all nucleons 
         noofnuc =  0
!               repeat until sum of mass becomes >= amn
!         *** until loop*** 
         do while (.true.)
            if( mno .ge. 10  .and. nooffragments .eq. 0) then
!                   for A >= 10, first sampled fragment must be
!                   non-nucleon to avoid too many nucleons.
                  call rndc(u)
                  u=(1.-cfragmentationtbl(ihg,2))*u +
     *                  cfragmentationtbl(ihg,2)
!                     find CfragmentationTbl >= u to sample fragrment
                  call kfrge(cfragmentationtbl(ihg,1), maxheavyg,
     *            ihg, u,  ihgf, jcon)
!                            sampled group index
            else
               call rndc(u)
!                    find FragmentaionTbl >= u to sample fragrment
               call kfrge(cfragmentationtbl(ihg,1), maxheavyg,
     *           ihg, u,  ihgf, jcon)
            endif
!
            if(ihgf .eq. 1)then
!               count nucleon
               noofnuc = noofnuc + 1
!               make ptcl;  nucleon charge is reset later
               call cmkptc(knuc, regptcl,
     *          1, nuc(noofnuc))
               nuc(noofnuc)%fm%p(4) = epn
            else
!               count fragment
               nooffragments=nooffragments+1
!               make ptcl; 
               call cmkptc(heavyg2code(ihgf), regptcl,
     *          1, fra(nooffragments))
!                 set energy
               fra(nooffragments)%fm%p(4) = heavyg2massn(ihgf) * epn
            endif
!               count mass no.
            msumf=msumf+heavyg2massn(ihgf)
            if  (msumf .ge. mn)
     *                      goto 50
         enddo
   50    continue
         if(msumf .le. mx) then
!           remaining ptcls should be nucleons
            do i = 1, mno - msumf
               noofnuc = noofnuc + 1
               call cmkptc(knuc, regptcl, 1, nuc(noofnuc) )
!                set energy
               nuc(noofnuc)%fm%p(4) = epn
            enddo
!             get sum of fragment charge
            zfrag = 0
            do i = 1, nooffragments
               zfrag = zfrag + fra(i)%charge
            enddo
            if(zfrag .le. pj%charge) goto 100
         endif
      enddo
  100 continue
!          reset nucleon charge
      call cresetnucchg(nuc, noofnuc, hchg - zfrag)
!
!               sample interacting nucleon no.
      call csampintenuc(pj, ia,  noofnuc, noofintenuc)

!           sample fragment mom. all frag. 
      call csampfragmom(fra,  nooffragments)
!           sample non-interacting nuc. mom.
      call csampfragmom(nuc, noofnuc - noofintenuc)
!           set interacting nuc. mom
      do   i=noofnuc, noofnuc - noofintenuc + 1, -1
         nuc(i)%fm%p(1) = 0.
         nuc(i)%fm%p(2) = 0.
         nuc(i)%fm%p(3) = sqrt(
     *      max( nuc(i)%fm%p(4)**2 - nuc(i)%mass**2, 0.d0) )
      enddo

Here is the call graph for this function:

csampfragmom()

subroutine csampfragmom	(	type(ptcl), dimension(nf)	a,
		integer	nf
	)

Definition at line 427 of file cheavyInt.f.

References csampfragpt(), d, kcossn(), and true.

Referenced by csampfragments().

      implicit none

#include  "Zptcl.h"
      integer nf
      type(ptcl):: a(nf)

      integer i, nc
      real*8 pt, p, cs, sn

       do   i=1, nf
!            sample fragment pt
         nc=0
         p=sqrt( a(i)%fm%p(4)**2- a(i)%mass**2 )
!         *** until loop*** 
         do while (.true.)
            call csampfragpt(a(i),  pt)
            nc=nc+1
            if         (pt .lt. p .or. nc .eq. 10)
     *                      goto 10
         enddo
   10    continue
         if(nc .ge. 10) then
             pt=min(1.d-10, p)
         endif
!               set pt and pz
         a(i)%fm%p(3) = sqrt(p**2-pt**2)
         call kcossn(cs, sn)
         a(i)%fm%p(1) = pt*cs
         a(i)%fm%p(2) = pt*sn
       enddo

Here is the call graph for this function:

Here is the caller graph for this function:

csampfragpt()

subroutine csampfragpt	(	type(ptcl)	aPtcl,
		real*8	pt
	)

Definition at line 469 of file cheavyInt.f.

References knuc, and rndc().

Referenced by csampfragmom().

      implicit none

#include  "Zptcl.h"
#include  "Zcode.h"
#include  "Zheavyp.h"
      real*8  pt      ! sampled pt in GeV/c

      type(ptcl):: aptcl
!
!
      real*8 avpt, u
!
      if(aptcl%code .eq. knuc) then
        avpt= ptavnonintenuc
      else
        avpt= ptavfrag
      endif
!
      call rndc(u)
!          note avpt is not <pt> but <pt>/sqrt(pi/2)
      pt  =  avpt *sqrt(- log(u)* 2 )

Here is the call graph for this function:

Here is the caller graph for this function:

csampintenuc()

subroutine csampintenuc	(	type(ptcl)	pj,
		integer	ia,
		integer	noOfNuc,
		integer	noOfInteNuc
	)

Definition at line 392 of file cheavyInt.f.

References caveintenuc(), and kbinom().

Referenced by csampfragments().

      implicit none

#include  "Zcode.h"
#include  "Zptcl.h"
#include  "Zheavyp.h"
      
      integer noofnuc, noofintenuc, ia  ! ia: # of target nuc. 
      type(ptcl):: pj  ! projectile heavy
      integer ihg
      real*8  avintn, avnn
!

!                heavy group index
      ihg=charge2heavyg(pj%charge)
!                get average no. of interacting nucleons
      call caveintenuc(pj, ia,  avintn)
!                average no. of nucleons in fragments
      avnn = fragmenttbl(ihg, 1)
!       sample interacting nucleon number by binormial distribution
!       with prob. avintn/avnn
      call kbinom( avintn/avnn, noofnuc, noofintenuc)

Here is the call graph for this function:

Here is the caller graph for this function: