cheavyInt.f

Go to the documentation of this file.

!     *****************************************************************
!     *                                                               *
!     *  cheavyInt: treats heavy particle interactions
!     *                                                               *
!     *****************************************************************
!
!  -- process --
!             1) samples fragment ptcls using fragmentation parameters
!                and determines the no. of interacting nucleons.
!             2) gives break-up angle for fragments other than interacting
!                nucleons
!             3) for interacting nucleons, calls  chAcol to make
!                multiple production.
!
        subroutine cheavyint(pj, ia, iz,  a, ntp)
#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
      end
!     *****************************************************************
!     *                                                               *
!     *  csampFragments: sample and set fragmentation ptcls           *
!     *                                                               *
!     *****************************************************************
!
!      /usage/  call csampFragments(pj, ia, fra, noOfFragments, 
!    *    nuc,  noOfNuc, noOfInteNuc)
!
!       1)  samples fragment ptcl one by one by referring cfptbl
!           until sum of fragment mass no. exceeds incident mass no.
!           - 3 so that no heavy ptcl can emerg more.  then the 
!           remaining ones
!           are assigned to nucleons, if any.  If the sum exceeds
!           incident mass no. during sampling, retrial is made from the
!           first.   For each sampled fragment, energy is given pro-
!           potionally to its mass.   charge, mass and kind are also
!           given according to cosmos convention.
!
!        2) fra is to store heavy fragments
!           nuc is to store all nucleons
!        3) charge of the nucleons is reset because process 1) assigns
!           only proton
!        4) samples no. of interacting nucleons
!        5) samples pt of fragments other than interacting nucleons
!           and convert it to ptx, pty
!
!      noOfFragments:  # of heavy fragments stored in fra
!      noOfNuc:       # of nucleons  stored in nuc.
!      noOfInteNuc:  # of interacting nucleons among them
!
!
        subroutine csampfragments(pj, ia,
     *       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
      end
!     *****************************************************************
!     *                                                               *
!     *  cresetNucChg:  reset charge of nucleons emerging from heavy *
!     *                                                               *
!     *****************************************************************
!
       subroutine cresetnucchg(nuc, noOfNuc, z)
!
       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
      end
!     *****************************************************************
!     *                                                               *
!     *  csampInteNuc:  sample interacting nucleon number                   *
!     *                                                               *
!     *****************************************************************
!
!        1) gets average no. of interacting nucleons
!        2) gets average no. of nucleons from heavy
!        3) using the ratio of 1)/2) and binomial distribution,
!           samples no. of interacting nucleons.

!
!                           =   =   =   =
!
      subroutine csampintenuc(pj, ia,  noOfNuc, noOfInteNuc)

      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)
      end
!     *****************************************************************
!     *                                                               *
!     *  csampFragMom:  sample  px, py, pz of fragments
!     *                                                               *
!     *****************************************************************
!
!      /usage/     call csampFragMom(a, nf)
!
!        pt is sampled by gaussian type distribution and stored in
!        'a'. ptx,pty are also stored.
!
      subroutine csampfragmom( a, nf )
      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
      end
!     *****************************************************************
!     *                                                               *
!     *  csampFragPt:   sample fragment pt                      
!     *                                                               *
!     *****************************************************************
!
!  -- process --
!        pt is sampled by gaussian type distribution:
!           exp(- x**2) dx**2  where x= 2pt/sqrt(pi)
!
      subroutine csampfragpt(aPtcl, pt)
      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 )
      end
!     *****************************************************************
!     *                                                               *
!     *  cmakeFragTbl:  makes fragmentation parameter table for 
!     *                 heavy ptcl frgmentation sampling              
!     *                                                               *
!     *****************************************************************
!
!        FragmentaionTbl(i,j) is assumed to have <no. of heavy of group 
!        j>  when heavy of group i fragments.  For each group i, 
!        FragmentTbl is normalized so that the sum of them
!        becomes 1 and then made to be cumulative table such that
!        CfragmentaionTbl(i,j) <= CragmentaionTbl(i,j+).
!
!

      subroutine cmakefragtbl
      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
      end
      subroutine chg2massn(hg, massnum)
      implicit none
#include "Zcode.h"
#include "Zheavyp.h"
      integer massnum, hg
      massnum = heavyg2massn(hg)
      end
      subroutine chg2charge(hg, charge)
      implicit none
#include "Zcode.h"
#include "Zheavyp.h"
      integer charge, hg
      charge = heavyg2charge(hg)
      end
      subroutine ccode2hvgrp(code, hg)
      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
      end
      subroutine ccode2mass(code, mass)
      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
      end