cgnlp.f File Reference

Go to the source code of this file.

Functions/Subroutines
subroutine	cgnlp (a, ntp, icon)
subroutine	ccpmul (roots, avn)
subroutine	ccylps (ntp, a, icon)
subroutine	c1pion (a, ntp, icon)
subroutine	caspti
subroutine	caspt (a, ntp)
subroutine	ciso (ntp, a, icon)
subroutine	cytoe (a, n)
subroutine	cspt (avpt, nptcl, a, ntpc)
subroutine	cptcns (a, nt, ptav)

Function/Subroutine Documentation

c1pion()

subroutine c1pion	(	type(ptcl), dimension(*)	a,
		integer	ntp,
		integer	icon
	)

Definition at line 142 of file cgnlp.f.

References cbst1(), cibst1(), cmkptc(), cnbdc3(), cnbdc4(), and kpion.

Referenced by cgnlp().

!         Method:
!          Now, the recoil leadings and missing mass and its
!          4 momentum are given.  We change this missing mass
!          to be a pion and make a transformaton so that
!          3 particles satify the 4 momentum conservation.            
!       a(1): /ptcl/. Outupt. produced pion
!        ntp: integer. Output.  1  if icon =0 (produced)
!                               0  if icon =1 (n.g)
!       icon: integer. Output.  0--> ok. 1--> n.g
!      ********************************************
       implicit none
!----       include '../../Zptcl.h'
#include  "Zptcl.h"
!----       include '../../Zmass.h'
#include  "Zmass.h"
!----       include '../Zevhnv.h'
#include  "Zevhnv.h"
#include  "Zcode.h"

!
       integer ntp, icon
       type(ptcl)::  a(*)

       type(ptcl):: p3(3)   ! to store 3 particles
       real*8 amu(3), roots, gzai
       integer  csum0, csum1, charge
       character*80  msg 
!
!       modify E,P following cpc p.364
!       get gzai for q~=gzai*p~
!       q0=sqrt(m**2 + gzai**2(p0**2-mu**2) )
!     
        p3(1) = rpjcms   ! reooil proj.  in cms
        p3(2) = rtgcms   ! recoil trgt.  in cms
        p3(3) = missingp ! missing mass in cms
        roots = cmsp%mass  

!           present mass
         amu(1)=pjlab%mass
         amu(2)=tglab%mass
         amu(3)=missingp%mass  ! missing mass
!            

         csum0 = pjlab%charge + tglab%charge
         csum1 = rpjcms%charge +rtgcms%charge
         if( abs(csum1-csum0) .gt. 1) then
!                retry once more
            icon =1 
            goto 900
         else
            charge = csum0 - csum1
!             true mass    ! we want missing --> pion mass.
!             modify missing mass to be the pion mass
            call cmkptc(kpion, 0, charge, p3(3))
         endif   
!             get convesion factor
!          note: in cnbdcy, amu is not used, because 
!                it is always 0. (start from 0 mass ptcl)
!                and transform to massive case (true mass
!                is in /ptcl/ data.)
!                However, here, the trial mass is missing mass 
!                and must be given in amu. True mass is
!                in /ptcl/ p3
         call cnbdc3(3, roots, p3, amu, 1,
     *   gzai, icon)
         if(icon .eq. 0) then
             call cnbdc4(3, p3,  amu, 1, gzai)
!            ________________________________
!            write(*,*) ' p3 sumpx',
!    *       (p3(1).fm.p(1)+p3(2).fm.p(1)+p3(3).fm.p(1))
!            write(*,*) ' p3 sumpy',
!    *       (p3(1).fm.p(2)+p3(2).fm.p(2)+p3(3).fm.p(2))
!            write(*,*) ' p3 sumpz',
!    *       (p3(1).fm.p(3)+p3(2).fm.p(3)+p3(3).fm.p(3))
!            write(*,*) ' p3 sum e (gev)',
!    *       ( p3(1).fm.p(4)+p3(2).fm.p(4)+p3(3).fm.p(4))
!            write(*,*) ' roots=',roots
!            ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!               reset leadings
             rpjcms = p3(1)
             rtgcms = p3(2)
!                boost  proj. recoil into lab (cms -> lab)
!                note:  tg(7).fm = cms 4 mom. in lab
             call cibst1(1, cmsp, rpjcms, rpjlab)
!                boost target recoil into lab
             call cibst1(2, cmsp, rtgcms, rtglab)
!               boost proj to target rest system. (may not be used)
             call cbst1(1, tglab, rpjlab, rpjtatr)
!               same for projectile rest system. (may not be used)
             call cbst1(1, pjlab, rtglab, rtgpatr)
!                move pion
             a(1) = p3(3)
             ntp = 1
         else
!            ____________________________________________
             write(0,*) ' failed to adjust  missing mass=',
     *       missingp%mass, ' into pion mass. '
!            ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
             ntp = 0
             icon = 1
         endif
 900     continue

Here is the call graph for this function:

Here is the caller graph for this function:

caspt()

subroutine caspt	(	type(ptcl), dimension(ntp)	a,
		integer	ntp
	)

Definition at line 292 of file cgnlp.f.

References cspt().

Referenced by ccylps().

!           generation of pt
!      *****************
!----        include '../../Zptcl.h'
#include  "Zptcl.h"
!----        include '../../Zcode.h'
#include  "Zcode.h"
!----        include '../Zevhnv.h'
#include  "Zevhnv.h"
        integer ntp
        type(ptcl):: a(ntp)

        real*8  pttmp
        integer ntpc
!           <>pt for NN~,          DD~
!        real*8 ptnnb/500.d-3/, ptddb/750.d-3/, ptavpi/330.d-3/,
!     *         ptavk/479.d-3/, ptaveta/480.d-3/
        real*8 ptnnb/250.d-3/, ptddb/370.d-3/, ptavpi/180.d-3/,
     *         ptavk/240.d-3/, ptaveta/240.d-3/
!        real*8 ptnnb/240.d-3/, ptddb/360.d-3/, ptavpi/170.d-3/,
!     *         ptavk/230.d-3/, ptaveta/230.d-3/
!
        ntpc = 0
!                 (n,n~)
        pttmp=ptnnb* ptnorm/ptavpi
!                     sample pt
        call cspt(pttmp, nnnb, a, ntpc)

!                (d,d~)
        pttmp = ptddb* ptnorm/ptavpi
        call cspt(pttmp, nddb, a, ntpc)
!                 pi+/-
        pttmp=ptnorm
        call cspt(pttmp, npic,  a, ntpc)
!                 pi0
        call cspt(pttmp, npi0,  a, ntpc)
!                eta
        pttmp = ptaveta*ptnorm/ptavpi
        call cspt(pttmp, neta, a, ntpc)
!                 kaon +/ -
        pttmp = ptavk*ptnorm/ptavpi
        call cspt(pttmp, nkch, a, ntpc)
!                 k0
        call cspt(pttmp, nk0,  a, ntpc)

Here is the call graph for this function:

Here is the caller graph for this function:

caspti()

subroutine caspti

(

)

Definition at line 247 of file cgnlp.f.

References d0, and masn.

Referenced by ccylps().

!     *********************************************************
!          preparation for assigning  pt
      implicit none
!----        include '../../Zptcl.h'
#include  "Zptcl.h"
!----        include '../../Zmass.h'
#include  "Zmass.h"
!----        include '../Zevhnv.h'
#include  "Zevhnv.h"
!
      real*8  efe0, xx
!      real*8  ptbase/330.d-3/
      real*8  ptbase/180.d-3/
!      real*8  ptbase/170.d-3/
!
!              some sort of effective E0, which I cannot
!           remember the reasoning now.
!        efe0 = (Efrs**2 + 2*(masn + Pjlab.mass)*Efrs)/2
!     *          /masn
        efe0 =max(1.d0,  (efrs**2 - masn*2- pjlab%mass**2)/2
     *          /masn)

!
!              Powerexp in   pt**pw2 * exp(-pt/..) dpt case.
        powerexp=1.0d0 + (efe0)**(-0.05d0)
        if(efe0 .lt. 40.) then
           xx = (efe0 - 40.)/5.5 + 6.
           ptnorm = ptbase *exp(xx)/(1+exp(xx))
        elseif(efe0 .lt. 40000.) then
           ptnorm= ptbase
        else
           ptnorm= ptbase
!           Ptnorm = ptbase *( efe0/40000.)**0.025
        endif
        probpower=min(0.04d0*nch, 0.33d0)  
        powerp=3.d0 +
     *        1.d0/
     *     (0.01 + 0.01*3.29* missingp%mass**0.3)
        if(powerp .gt. 100.d0) then
             probpower=0.d0
        endif

Here is the caller graph for this function:

ccpmul()

subroutine ccpmul	(	real*8	roots,
		real*8	avn
	)

Definition at line 69 of file cgnlp.f.

Referenced by csnchp().

!     *****************************
!      average charged particle multiplicity at root(s)
!     /**** UA5 parameterization ***/
!      root(s)  GeV
       implicit none
#include "Zevhnp.h"

       real*8 roots, avn
!       real*8 lambda/0.3/, no/0.6135/, qcdErg/1000./  ! up to 1000 GeV.
!                              | this is wrong (total  multiplicity)
        real*8 lambda/0.3/, no/0.34/, qcderg/1000./  
!                           
!               e UA5 data.
!                   original
!       avn= 7.2* roots**(2*0.127) -7.   including  leading.
!
         avn = (no* exp(sqrt(23./18. * log( (roots/lambda))*2))
     *          -3.5)*0.8 
         if(avn .lt. 0.1) then
            avn = 0.1
         endif
!
!

Here is the caller graph for this function:

ccylps()

subroutine ccylps	(	integer	ntp,
		type(ptcl), dimension(ntp)	a,
		integer	icon
	)

Definition at line 99 of file cgnlp.f.

References caspt(), caspti(), cgrap(), cibst1(), cptcns(), csptxy(), and cytoe().

Referenced by cgnlp().

       implicit none
!----       include  '../../Zptcl.h'
#include  "Zptcl.h"
!----       include  '../Zevhnv.h'
#include  "Zevhnv.h"
       integer ntp, icon
       type(ptcl)::  a(ntp)
!
       real*8 ptav, w
       integer i
!
!          assign pt 
!              init for const.
         call caspti
!                            loc pt = loc pz
         call caspt(a, ntp)
!          pt---> ptx,pty
         call csptxy(a, ntp)
!            forced conservation of pt.
         call cptcns(a, ntp, ptav)
!           generation of rapidity for missing mass
         w = missingp%mass
!                             loc tm = loc pz
!                             loc rap = loc e
         call cgrap(w, ptav, ntp,  a, icon)
         if(icon .eq. 0) then
!             convert y into to missing mass system energy
            call cytoe(a, ntp)
!           ________________________________
!           call cchk(' after cytoe', a, ntp)
!           ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!                boosting ntp ptcls to cms
             do   i=1, ntp
                call cibst1(i, missingp, a(i), a(i))
             enddo
         endif

Here is the call graph for this function:

Here is the caller graph for this function:

cgnlp()

subroutine cgnlp	(	type(ptcl), dimension(*)	a,
		integer	ntp,
		integer	icon
	)

Definition at line 4 of file cgnlp.f.

References c1pion(), ccylps(), cfnptc(), ciso(), csnchp(), and true.

Referenced by chncol2().

!     ********************************************************
!            a: /ptcl/    Output. to get produced particles in cms
!          ntp:  output.  Integer. The number of particle in a.
!         icon:  Integer. Output.0--> o.k
!                                1--> n.g
       implicit none
#include  "Zptcl.h"
#include  "Zmass.h"
#include  "Zevhnv.h"
!
       type(ptcl):: a(*)
       integer ntp, icon
!      -------------------------------- 
       real*8   missgm, roots
       integer jcon
!


!       integer maxiso/10/  ! max # of ptcl from isotropic p.s   
       integer maxiso/0/  ! max # of ptcl from isotropic p%s   
!       logical ok
!
!              sample # of charged ptcls
!
       missgm = missingp%mass
       roots =  cmsp%mass
!        *** until loop*** 

       do while (.true.)
!                get average # and sampled # of charged ptcl
!                Avncharged & Nch
          call csnchp(jcon)
          if(jcon .ne. 0) then
             icon=1
             goto 900           ! return
          endif
!              fix # of pi+-, pi0, k+-,k0, nn~, dd~,
!              and give ptcl mass and code in a;
!              # of pi+- etc is obtained by calling
!              cqnptc(code, charge, nout)
          call cfnptc(a,  ntp)
           if(ntp .ge. 1 ) goto 10
!           if(ntp .eq. 1 .and. missgm .lt. 15.*maspic) goto 10
!           if(ntp .eq. 2 .and. missgm .lt. 30.*maspic) goto 10
       enddo
 10    continue
       if(ntp .le. maxiso .and. ntp .ge. 2) then
!               use isotripc p.s
            call ciso(ntp, a, icon)
            if(icon .eq. 1) then
!                   <Pt> is too large so use cylindrical p.s
               call ccylps(ntp, a, icon)
            elseif(icon .eq. 2) then
                icon=1
            endif
        elseif(ntp .ge. 2) then
           call ccylps(ntp,  a, icon)
        else
!              assume 1 (or 0) pion production in cms
            call c1pion(a, ntp, icon)
        endif
  900  continue

Here is the call graph for this function:

Here is the caller graph for this function:

ciso()

subroutine ciso	(	integer	ntp,
		type(ptcl), dimension(ntp)	a,
		integer	icon
	)

Definition at line 338 of file cgnlp.f.

References cerrormsg(), cibst1(), cnbdcy(), d0, and true.

Referenced by cgnlp().

!           isotropic phase space
!          icon: 0---> o.k
!                1---> <pt> is big,   better to use cylindrical p.s
!                2---> n.g
!      ****************************************************
       implicit none
!----       include  '../../Zptcl.h'
#include  "Zptcl.h"
!----       include  '../Zevhnv.h'
#include  "Zevhnv.h"
       integer ntp, icon
       type(ptcl)::   a(ntp)
!
       integer nc, jcon, i
       real*8 ptc1/500.d-3/, ptc2/700.d-3/
!
       real*8 ptlm1, ptlm2, missgm, wg, sumpt, avpt
!
       character*80  msg

       missgm = missingp%mass
       nc=0
       if(pjlab%fm%p(4) .lt. 1000.d0) then
             ptlm1=ptc1*(pjlab%fm%p(4)/1000.d0)**0.04
             ptlm2=ptc2*(pjlab%fm%p(4)/1000.d0)**0.04
       else
             ptlm1=ptc1
             ptlm2=ptc2
       endif
       if(3.1415/4.0* missgm/ntp .gt. ptlm1) then
             jcon=1
       else
!             *** until loop*** 
             do while (.true.)
!                  isotropic p.s decay; take almost all weights
!                  wg=.95 does not make any diff.(slower speed)
!                call cnbdcy(ntp, missgm, a, 1, wg, jcon)
                call cnbdcy(ntp, missgm, a, 0, wg, jcon)
                nc=nc+1
                if(nc .gt. 15) then
                    jcon=3
                    write(msg,*) ' cnbdcy fail but try further'
                    call cerrormsg(msg, 1)
                endif
                if(jcon .ne. 0 .or. wg .gt. .05d0) goto 99
             enddo
   99        continue
             if(jcon .eq. 0) then
!                 boost 4 momentum into cms
                 do i=1, ntp
                    call cibst1(i, missingp, a(i), a(i))
                 enddo
!                   get <pt> and see if it is too large
!                        get sum of  pt
                 sumpt=0.d0
                 do   i=1, ntp
                     sumpt=sumpt +
     *                  sqrt(a(i)%fm%p(1)**2+a(i)%fm%p(2)**2)
                 enddo
                 avpt=sumpt/ntp
                 if(avpt .gt. ptlm2) then
                    jcon = 2
                 elseif(avpt .gt. ptlm1) then
                    jcon = 1
                 else
                    jcon = 0
                 endif
             else
                 jcon = 2
             endif
         endif
         icon=jcon

Here is the call graph for this function:

Here is the caller graph for this function:

cptcns()

subroutine cptcns	(	type(ptcl), dimension(nt)	a,
		integer	nt,
		real*8	ptav
	)

Definition at line 492 of file cgnlp.f.

References d, and d0.

Referenced by ccylps().

       implicit none
!----       include '../../Zptcl.h'
#include  "Zptcl.h"
       integer nt
       type(ptcl):: a(nt)
       real*8 ptav 
!
       real*8 sumpx, sumpy, sumpt, cfx, cfy
       real*8 sumpt2, cf
       integer i 
!
      if(nt .ge. 2) then 
!           get sum of pt, ptx, pty
         sumpx = 0.d0
         sumpy = 0.d0
         sumpt = 0.d0
         do i=1, nt
            sumpt = sumpt + a(i)%fm%p(3)
            sumpx = sumpx + a(i)%fm%p(1)
            sumpy = sumpy + a(i)%fm%p(2)
         enddo
         if(sumpt .gt. 0.d0) then 
!              correction factor
            cfx=sumpx/sumpt
            cfy=sumpy/sumpt 
!
            sumpt2=0.d0
            do i=1, nt
               a(i)%fm%p(1) = a(i)%fm%p(1) - a(i)%fm%p(3) * cfx
               a(i)%fm%p(2) = a(i)%fm%p(2) - a(i)%fm%p(3) * cfy
               a(i)%fm%p(3) = sqrt(a(i)%fm%p(1)**2 + a(i)%fm%p(2)**2 )
               sumpt2=a(i)%fm%p(3) + sumpt2
            enddo 
!                adjust: sum (pt) = original value
            cf = sumpt/sumpt2 
!                  multipliy cf to pt and ptx,pty
            do i=1, nt
                a(i)%fm%p(3) = a(i)%fm%p(3)*cf
                a(i)%fm%p(1) = a(i)%fm%p(1)*cf
                a(i)%fm%p(2) = a(i)%fm%p(2)*cf
            enddo
            ptav=sumpt/nt
        else
            ptav = 1.d-1
        endif
      elseif(nt .gt. 0) then
         ptav=a(1)%fm%p(3)
      else
         ptav = 1.d-1
      endif

Here is the caller graph for this function:

cspt()

subroutine cspt	(	real*8	avpt,
		integer	nptcl,
		type(ptcl), dimension(*)	a,
		integer	ntpc
	)

Definition at line 434 of file cgnlp.f.

References cspwpt(), ksgmrm(), and rndc().

Referenced by caspt().

! to sample a pt 
! avpt: real*8. Input. average pt in exponential part.  
! nptcl: integer. Input. # of ptcls to be assigned
!     a: /ptcl/.  output. a.fm.p(3) is given a pt value
!  ntpc: integer. input/oututp. a(ntpc+1) is the first
!                 ptcl pos.  ntpc is incremented by nptcl
!                 on return.

        implicit none
!----        include '../../Zptcl.h'
#include  "Zptcl.h"
!----        include '../Zevhnv.h'
#include  "Zevhnv.h"
        integer nptcl, ntpc
        type(ptcl):: a(*)
        real*8 avpt
!
!        real*8 u, bpt/1.7d0/, pt
!        real*8 u, bpt/1.5d0/, pt
        real*8 u, bpt/1.0d0/, pt


        integer nc
        
        do nc = 1, nptcl
           call rndc(u)
           if(u .lt. probpower) then 
!                  power pt
                 call cspwpt(bpt, powerp, pt)
           else 
!                   pt**pw2 * exp(-pt/..)dpt

                call ksgmrm(powerexp, avpt, pt)
           endif
           ntpc = ntpc +1
           a(ntpc)%fm%p(3) = pt
        enddo

Here is the call graph for this function:

Here is the caller graph for this function:

cytoe()

subroutine cytoe	(	type(ptcl), dimension(n)	a,
		integer	n
	)

Definition at line 415 of file cgnlp.f.

Referenced by ccylps().

      implicit none
!----       include  '../../Zptcl.h'
#include  "Zptcl.h"
       integer n
       type(ptcl):: a(n)
!
       integer i
       real*8 etemp
          do   i=1, n
!              note that rap and e have the same pos.
!             etemp=a(i).fm.tm *  cosh( a(i).fm.rap  )
             etemp=a(i)%fm%p(3)*  cosh( a(i)%fm%p(4)  )
!             a(i).fm.p(3) = a(i).fm.tm * sinh( a(i).fm.rap )
             a(i)%fm%p(3) = a(i)%fm%p(3)* sinh( a(i)%fm%p(4))
             a(i)%fm%p(4) = etemp
          enddo

Here is the caller graph for this function: