#include "Zintmodel.h"

#ifdef MODIFYX
#include "csoftenPiK.f"      
#endif

#include "ZcosmosBD.h"
#ifdef MODIFYX
      use SoftenPiK
#endif
      implicit none
#include  "Zptcl.h"
#include  "Ztrackp.h"
#include  "Zevhnp.h"

      include "Zprivate.h"

      integer i, nev, j, ntp, eof, ntpo
      record /ptcl/ w(maxn)

      call init
      do j = 1, nevent
         if(inpfileno .gt. 0) then
            call readinpfile(eof)
            if(eof .eq. 1) then
               write(0,*) 
     *        ' number of events generated is ',j-1
               goto 100
            endif
            call formpjtg(0)
         endif
         call gencol(w, ntp)
#ifdef  MODIFYX
c                 in cms we work
         call csoftenFE(pj, fwbw,  w, ntp, ntpo)
         ntp = ntpo
#endif
         call cutbyangle(w, ntp, ntpo)
         ntp = ntpo
         call sortbyke(w, ntp)  ! sort by kinetic energy 
         if(Trace .gt. 0) then
            call outtrace(j, w, ntp)
         endif
         if(outresul1 == 0 ) then
            call outresul(w, ntp)
         elseif(outresul1 == 1 ) then
            call outresulB(w, ntp)
         endif
      enddo
      write(0,*) 
     *  ' number of events generated is ',nevent
 100  continue
      write(0,*)
     * "Equivalent Lab. Energy was ", LabEquivE," GeV/n"
      end

      subroutine init
      implicit none
#include  "Zptcl.h"
#include  "Zcode.h"
#include  "Zevhnp.h"
#include  "Zevhnv.h"
#include  "Zmass.h"
#include  "Zmanager.h"
#include  "Zmanagerp.h"
#include  "Ztrackp.h"

      include  "Zprivate.h"
      character*200 input, file
      character*20 uid
      integer klena, icon, eof

      external cblkManager
      external cblkEvhnp

      call creadParam(5)
      

      if(TraceDir .eq. ' ') then
         call cgetLoginN(uid)
         TraceDir = '/tmp/'//uid(1:klena(uid))
      endif

      if(DestEventNo(2) .eq. 0) then
         nevent =abs( DestEventNo(1) )
      else
         nevent = abs( DestEventNo(2) )
      endif
      call cmkSeed(InitRN(1), InitRN)
      call rnd1i(InitRN)        ! random number init.
      call cqUhookr(1, wzmin)
      call cqUhookr(2, wzmax)
      write(0,*) ' cos* min  max=', wzmin, wzmax
      call cqUhookr(3, trackl)
      call cqUhooki(1, outzero)
      call cqUhooki(2, outresul1)
c       make projectile going +z
      call cqUhookc(1, input) 
      if(input(1:4) .eq. "file") then
         read(input(5:10), *) inpfileno
         xyz=index(input, "xyz")
         call cqUhookc(2, input) 
         file = ' '
         file=input(1:klena(input))
         call copenfw2(inpfileno, file, 1, icon)
         if(icon .ne. 1) then
            write(0,*)
     *      ' file=', file, ' cannot be opened in Gencol'
            stop 9999
         endif
         call readinpfile(eof)
c          once rewind to read successively for event generation
         rewind inpfileno
      else
         inpfileno=0
         read(input, *) 
     *        pjcode, pjsub, pjchg, pjpx, pjpy, pjpz  !  proj. mom/n
         call cqUhookc(2, input) 
         input = trim(input)//" /"
         read(input, *) 
     *        tgcode, tgsub, tgchg, tgpx, tgpy, tgpz, targetName  ! target. mom/n
         call cqUhookc(3, input)
         if(input .ne. ' ') then
c            read(input, *) xpos, ypos, zpos
            read(input, *) xpos(1:3)
            xyz = 1
         else
            xyz = 0
         endif
      endif
      call formpjtg(1)    ! form proj. and target
      call cfixPrefix('configDummy')
      call csetCosOrEpi('gencol') ! this is probably not used


      if( index( IntModel,'qgsjet1') .ne. 0 ) then
#ifdef QGSJET1
         call qgs01init
         ActiveMdl = 'qgsjet1'
#else
         write(0,*) 'to use qgsjet1,  define it  in Zintmodel.h'
#endif
      elseif(index (IntModel, 'sibyll') .ne. 0 )  then
ccc  #ifdef  SIBYLL
         call csibyllinit
         ActiveMdl = 'sibyll'
cc #else
ccc        write(0,*) 'to use sibyll, define it in Zintmodel.h'
ccc #endif
      else
         call cintModels('gencol')
         call cfixModel( plab )
      endif

      write(0, *) 'Active int. model=',ActiveMdl

      if( plab.code == kgnuc ) then
         LabEquivE= plab.fm.p(4)/plab.subcode
      else
         LabEquivE=plab.fm.p(4)
      endif
      if( ActiveMdl == "dpmjet3" ) then
         call checkEnergy
         call checkPtcls
      endif
!             for muon N.I, we must read  mutab for which we must read 
!         all  elemag table too.
      if(pj.code == kmuon ) then
cc         call epReadMediaForMuon  !////////////
      endif
      end
      subroutine  epReadMediaForMuon
#include  "Zcode.h"
#include  "Zptcl.h"
      
      include  "Zprivate.h"
      integer icon

      if( targetName /= " " ) then
         MediaNo = 1
cccc         call eprdMFile(targetName, icon)  !!////////////
         if( icon /= 0 ) then
            write(0,*) " targetName =", targetName
            write(0,*) " not acceptable "
            stop
         endif
      else
         write(0,*) " projectile code =", plab.code, "= muon"
         write(0,*) "   In this case, media name of the target must"
         write(0,*) " follow the target partcile momentum, e.g,"
         write(0,*) " UserHookc = '1  -1 0   0. 0     500',"
         write(0,*) "             '9  16 8   0  0   0. LiqO2',"
         write(0,*) "             "
         stop
      endif                       
      end


      subroutine checkEnergy
      implicit none
#include "Zptcl.h"
#include "Zprivate.h"
#include "Zcode.h"
      integer icon

      character(80):: input
      character(8):: bb(3)
      integer:: nr
      real(8):: Einput

      if( pj.code == kphoton ) then
         if( tg.code == kgnuc ) then
            if( pj.fm.p(4) < 7. ) then
               write(0,*)
     *         'photon E must be >= 7 for A target'
               stop
            endif
         else
            if( pj.fm.p(4) < 6. ) then
               write(0,*)
     *         'photon E must be >= 6 for non-A target'
               stop
            endif
         endif
      endif
         
      call copenf(11, "dpmjet.inp", icon)
      if(icon /= 0 ) then
         write(0,*) ' dpmjet.inp  cannot be opened'
         stop
      endif
      do while(.true.)
         read(11,'(a)') input
         call ksplit(input,8, 3, bb, nr)
         if( bb(1) == "ENERGY") then
            read(bb(2),*) Einput
            write(0,*) "ENERGY in dpmjet.inp =",  Einput
            if( Einput < LabEquivE) then
               write(0,*) " is  too small; give a value" 
               write(0,*) "little bit larger than ",LabEquivE
               write(0,*) "   if taget is A, ~1.5 % "
               write(0,*) "   if taget is p/n, very little"
                  write(0,*) 
     *           "E.g ",LabEquivE*1.001, " with ~ 4 digit accuracy" 
                  stop
!            elseif( Einput > LabEquivE*1.002 ) then
            elseif(  tg.code == kgnuc .and. 
     *              Einput > LabEquivE*1.040 ) then
               write(0,*) " is  too large; give a vlaue" 
               write(0,*) "little bit larger than ",LabEquivE
               write(0,*) 
     *           "E.g ",LabEquivE*1.03, " with ~ 4 digit accuracy" 
               stop
            elseif( tg.code /= kgnuc .and.
     *             Einput > LabEquivE*1.002 )  then
               write(0,*) " is  too large; give a vlaue" 
               write(0,*) "little bit larger than ",LabEquivE
               write(0,*) 
     *           "E.g ",LabEquivE*1.001, " with ~ 4 digit accuracy" 
               stop
            else
               write(0,*) " is close to ", LabEquivE, " and OK"
            endif
            close(11)
            exit
         endif
      enddo
      end

      subroutine checkPtcls
!       check projectile and target  ;
!       Are they same in param and dpmjet.inp ? 
      implicit none
#include "Zptcl.h"
#include "Zprivate.h"
#include "Zcode.h"
      integer icon

      character(80):: input
      character(8):: bb(5)
      integer:: nr
      character(8)::btype
      integer code, subc, charge

      call copenf(11, "dpmjet.inp", icon)
      if(icon /= 0 ) then
         write(0,*) ' dpmjet.inp  cannot be opened'
         stop
      endif
      do while(.true.)
         read(11,'(a)', end=100) input
         call ksplit(input,8, 5, bb, nr)

         if( bb(1) == "PROJPAR") then
            if( nr == 2 ) then
               btype=bb(2)
               call cbtype2cos(btype, code, subc, charge)
               if( code /= pj.code  .or. charge /= pj.charge ) then
                     ! subcode not checked 
                  write(0,*) 
     *                 ' projectile mismatch bw dpmjet.inp and param'
                  write(0,*) ' dpmjet.inp :', btype
                  write(0,*) ' param: ',  pj.code, pj.charge
                  stop
               endif
            elseif( nr == 3 ) then
               read( bb(2), *) subc
               read( bb(3), *) charge
               if( subc > 1 ) then
                  if( subc /= pj.subcode .or. charge /= pj.charge ) then
                     write(0,*) 'proj  mismatch bw dpmjet.inp and param'
                     stop
                  endif
               else
                  if( pj.code /= knuc .and. charge /= pj.charge ) then
                     write(0,*) 'proj  mismatch bw dpmjet.inp and param'
                     stop
                  endif
               endif
            endif
         elseif(bb(1) == "TARGETPAR") then
            if( nr == 2 ) then
               btype=bb(2)
               call cbtype2cos(btype, code, subc, charge)
               if( code /= tg.code  .or. charge /= tg.charge ) then
                     ! subcode not checked 
                  write(0,*) 
     *         ' target mismatch bw dpmjet.inp and param'
                  write(0,*) ' dpmjet.inp :', btype
                  write(0,*) ' param: ',  tg.code, tg.charge
                  stop
               endif
            elseif( nr == 3 ) then
               read( bb(2), *) subc
               read( bb(3), *) charge
               if( subc > 1 ) then
                  if( subc /= tg.subcode .or. charge /= tg.charge ) then
                     write(0,*)
     *                  'target mismatch bw dpmjet.inp and param'
                     stop
                  endif
               else
                  if( tg.code /= knuc .and. charge /= tg.charge ) then
                     write(0,*) 
     *                  'target mismatch bw dpmjet.inp and param'
                     stop
                  endif
               endif
            endif 
         endif
      enddo
 100  continue
      close(11)
      write(0,*)
     *  ' proj targ consistency btw dpmjet.inp and param: OK'
      end
      subroutine  cbtype2cos(btype, code, subc, charge)
      implicit none
#include "Zcode.h"
      character(8),intent(in):: btype
      integer,intent(out):: code, subc, charge
      
      select case(btype)
         case("PROTON") 
           code = knuc
           charge = 1
           subc = -1
        case("APROTON") 
           code = knuc
           charge = -1
           subc = 1
        case("PHOTON")
           code = kphoton
           charge = 0
           subc = 0
        case("NEUTRON")
           code = knuc
           charge = 0
           subc = -1
        case("ANEUTRON") 
           code = knuc
           charge = 0
           subc = 1
        case("PION+")
           code = kpion
           charge = 1
           subc = -1
        case("PION-")
           code = kpion
           charge = -1
           subc = 1
        case("KAON+")
           code = kkaon
           charge = 1
           subc = -1
        case("KAON-")
           code = kkaon
           charge = -1
           subc = 1
        case("LAMBDA")
           code = klambda
           charge = 0
           subc = -1
        case("ALAMBDA")
           code = klambda
           charge = 0
           subc = 1
        case("PIZERO")
           code = kpion
           charge = 0
           subc = 0
        case default
           write(0,*) ' not acceptable ptcl=',btype
           stop
        end select
!     'KAONSHRT' ,
!     &'SIGMA-  ' , 'SIGMA+  ' , 'SIGMAZER' , 'PIZERO  ' , 'KAONZERO' ,
!     &'AKAONZER' 
      end subroutine cbtype2cos

      subroutine readinpfile(eof)
      implicit none
#include "Zptcl.h"      
      include "Zprivate.h"

      integer eof ! output . data read--> 0
                  !   eof reached --> 1
      read(inpfileno,*, end=100)
     *     pjcode, pjsub, pjchg, pjpx, pjpy, pjpz
      read(inpfileno,*, end=100)
     *     tgcode, tgsub, tgchg, tgpx, tgpy, tgpz
      if(xyz .gt. 0 ) then
c         read(inpfileno,*, end=100) xpos, ypos, zpos
         read(inpfileno,*, end=100) xpos(1:3)
      endif
      eof = 0
      return
 100  continue
      eof = 1
      end
c     *******************
      subroutine formpjtg(confirm)
c     ******************
      implicit none
#include  "Zptcl.h"
#include  "Zcode.h"
#include  "Zevhnp.h"
#include  "Zevhnv.h"
#include  "Zmass.h"
#include  "Zmanager.h"
#include  "Zmanagerp.h"
#include  "Ztrackp.h"

      include "Zprivate.h"
      record /ptcl/resttg, Zcmsp, Zplab
      integer confirm  !  input. if 0, root s is  not printed.
                      !         else  printed
      real*8   roots, s
c         form projectile  and target
      integer::icon

      call cmkptc(pjcode, pjsub, pjchg,   pj)
c       pjmnum:    proj. mass number in integer
      if(pjcode .ne. kgnuc) then
         pjmnum = 1
      else
         pjmnum = pjsub
      endif
      pj.fm.p(1) = pjpx*pjmnum     ! total mom.
      pj.fm.p(2) = pjpy*pjmnum
      pj.fm.p(3) = pjpz*pjmnum
      pj.fm.p(4) =
     *     sqrt(pj.fm.p(1)**2 + pj.fm.p(2)**2 + pj.fm.p(3)**2
     *     + pj.mass**2)
         
c       make taget (rest or moving -z or ...)
      call cmkptc(tgcode, tgsub, tgchg, tg)
      if(tgcode .ne. kgnuc) then
         tgmnum = 1
      else
         tgmnum = tgsub
      endif
      tg.fm.p(1) = tgpx*tgmnum   ! total mom.
      tg.fm.p(2) = tgpy*tgmnum
      tg.fm.p(3) = tgpz*tgmnum
      tg.fm.p(4) =
     *     sqrt(tg.fm.p(1)**2 + tg.fm.p(2)**2 + tg.fm.p(3)**2 
     *   +   tg.mass**2)
c       
      if(tgpx .eq. 0. .and. tgpy .eq. 0. .and.
     *        tgpz .eq. 0.)  then
c     target is at rest;  s = (Ep+Et)^2 - (Pp+Pt)^2
c                           = (Ep+Mt)^2 - Pp^2
c                           =  Mp^2 +2EpMt +Mt^2
c   
         s= 2*pj.fm.p(4)*tg.mass +tg.mass**2 + pj.mass**2
      else
c         by  general formula;
c               Mp^2 + Mt^2 +2(Ep*Et - Pp.Pt)
         s = pj.mass**2 + tg.mass**2 +
     *       2*(pj.fm.p(4)*tg.fm.p(4) -
     *         dot_product(pj.fm.p(1:3), tg.fm.p(1:3) ) )
      endif
      roots = sqrt(s)
c      if(confirm .ne. 0) then
         write(0, *) ' roots/2=', roots/2
         write(0, *) 's,roots above are total not by /n'
c      endif
cc           boost to target rest system
      call cbst1(1, tg, pj, plab)
c         make rest target
      call cbst1(2, tg, tg, resttg)
       !   get /n cms
      call cgeqm2(plab, resttg, Cmsp, icon)   
c////////////
      write(0,*) ' Cmsp= (/n)',Cmsp.fm.p(:), Cmsp.mass
      write(0,*) ' plab=',plab.fm.p(4)
      write(0,*) ' tg=', tg.fm.p(4)
c///////////
      if(icon /= 0 ) then
         write(0,*) ' input 2 ptcls cannot form CMS'
         stop
      endif

!       Next part  will be used only by EPOS 
!       get equiv. CMS particle (Zcmsp) when projectile Pt is 0.
!       If crossing angle is /= 0, Cmsp and Zcmsp  
!       differ and we must use Zcmsp here!!!!.   2016/01/24
      zplab = plab
      zplab.fm.p(1:2) = 0.
      zplab.fm.p(3) = sqrt( zplab.fm.p(4)**2- zplab.mass**2)
       !   get /n cms
      call cgeqm2(zplab, resttg, Zcmsp, icon)   
c         inform Zcmsp to cintModels.f 
!      call cputGencolCMS(Cmsp) 
      call cputGencolCMS(Zcmsp) 
      end
c     ************************
      subroutine outresul(a, ntp)
      implicit none
#include  "Zptcl.h"
#include  "Zcode.h"
#include  "Zevhnp.h"
#include  "Zevhnv.h"
#include  "Zmass.h"
#include  "Zmanagerp.h"
#include  "Zmanager.h"
#include  "Ztrackp.h"
      include  "Zprivate.h"

      integer ntp
      record /ptcl/ a(ntp)
      integer  i, j
      real*8  p
      real*8  wx(3)
      integer nw, difcode(20)

      call getDiffCode(nw, difcode)

      do j = 1, ntp
         i = indx(j)
c         p= sqrt( a(i).fm.p(1)**2 + a(i).fm.p(2)**2
c     *        +      a(i).fm.p(3)**2 )
         p = sqrt(
     *  dot_product( a(i).fm.p(1:3), a(i).fm.p(1:3)) )
         wx(1:3)= a(i).fm.p(1:3)/p
c         wx = a(i).fm.p(1)/p               
c         wy = a(i).fm.p(2)/p               
c         wz = a(i).fm.p(3)/p               
         if(xyz .eq. 0) then
            write(*,
     *     '(i3,i5,i4,1p, g14.5,0p, 3f17.13, i8)')
     *        a(i).code, a(i).subcode, a(i).charge,
     *        a(i).fm.p(4)-a(i).mass, wx(1:3), j

!            write(*,
!     *     '(i3,i5,i4,1p, g14.5,0p, 3f17.13, 1p,3g12.4, i8)')
!     *        a(i).code, a(i).subcode, a(i).charge,
!     *        a(i).fm.p(4)-a(i).mass, wx(1:3), 
!     *        a(i).fm.p(1:3),   j

         else
            write(*,'(i3,i5,i4,g14.5,1p3E11.3,0p3f17.13, i8)')
     *        a(i).code, a(i).subcode, a(i).charge,
c     *        a(i).fm.p(4)-a(i).mass, xpos, ypos, zpos, 
     *        a(i).fm.p(4)-a(i).mass, xpos(1:3), wx(1:3), j
c     *        wx, wy, wz, j
         endif
      enddo
      if(ntp .gt. 0 .or. outzero .eq. 1) then
         write(*, *) 
      endif
      end

      subroutine  gencol(a, ntp)
      implicit none
#include  "Zptcl.h"
#include  "Zcode.h"
#include  "Zevhnv.h"
#include  "Zevhnp.h"
#include  "Zmanagerp.h"
      include "Zprivate.h"

      record /ptcl/  a(*)
c             projectile and target information (both befor
c             and after collision ) in different system.
c
      integer  ntp
      integer j
      integer tZ, tA
      real*8  xs
c     
      tA = tgmnum
      tZ =  tg.charge
      if(pj.code == kmuon ) then
!         sepcial treatment
         call epGencolByMuon(plab, tA, tZ,  a, ntp)
!!      elseif( pj.code == kphoton ) then
!         sepcial treatment
!!         call epGencolByPhoton(plab, tA, tZ,  a, ntp)
      elseif(pj.code == kelec .or.
     *      ( pj.code /= kgnuc .and. pj.code >knuc ) ) then
         write(0,*) ' ptcl code =', pj.code, ' not '
         write(0,*) ' supported in Gencol'
         stop
      else
!           init for 1 event 
         if(ActiveMdl .eq. 'qgsjet2' ) then
            call cxsecQGS(plab, tA, xs)
         elseif(ActiveMdl .eq. 'epos' ) then
            call ceposIniOneEvent(plab, tg, xs)
         endif

!           generate event
         if(ActiveMdl .eq. 'qgsjet1') then
#ifdef QGSJET1
            call qgs01event(plab, tA, tZ, a, ntp)
#endif
         elseif(ActiveMdl .eq. 'sibyll') then
ccc #ifdef SIBYLL
            call csibyllevent(plab, tA, tZ, a, ntp)
ccc #endif
         else
            call chAcol(plab, tA, tZ,  a, ntp)
         endif
      endif
      call getImpactParam(b)
      do j = 1, ntp
c               boost to original target system
         call cibst1(j, tg,  a(j), a(j))
      enddo
      end
      subroutine  epGencolByMuon
      implicit none
!      call epInfoPhotoP(incGp);  incGpをHowPhotopに焼き直し
!      call ep2cosCond
!      call cfixModel( cTrack.p )
!      call ciniSmpIntL   ! not related.
!      call epsmpNEPIntL(Media(MediaNo))  !st decay length only
!      call ep2cosCondr  !set FromEpics = .false.
!  if(   cTrack.p.fm.p(4) .gt.  Media(MediaNo).cnst.muNEmina )
!      call epmuNsmpP( Media(MediaNo),
!     *    cTrack.p.fm.p(4), prob, path)
!      epfixProc    !  ProcessNo 
!
!      epmuinte:
!       call epmuNsmpE(Media(MediaNo), Move.Track.p.fm.p(4),
!     *                   Et)
!      cTrack.p.fm.p(4) = cTrack.p.fm.p(4) - Et -->muon
!        if(Et .gt. 152.d-3) then
!            if(Media(MediaNo).mu.MuNI .eq. 3 ) then
!c             generate gamma-N interaction; employ gamma interaction  
!c             routine                                                 
!               cTrack.p.fm.p(4) = Et
!               call cmkptc( kphoton, 0, 0, cTrack.p)
!               call epe2p(cTrack) ! adjust momentum                   
!c                                                                     
!               call ep2cosPtcl( cTrack.p )
!c                 for small basic cross section case.                 
!               call epfixTarget2(ActiveMdl, Media(MediaNo))
!               call ep2cosCond2(Media(MediaNo).colA,
!     *                          Media(MediaNo).colZ)
!               call cphotop     ! Cosmos function     
!         call eppushPtcl(cTrack) ! use pos. info from this ptcl

      end

      subroutine  epGencolByPhoton(plab, tA, tZ,  a, ntp)
      implicit none
c  #include  "Zcode.h"
#include  "Ztrack.h"
#include  "Ztrackv.h"
#include  "Ztrackp.h"
#include  "Zelemagp.h"

c       include "Zprivate.h"

      record /ptcl/ plab  ! input proj. in Lab
      integer,intent(in):: tA   ! target A
      integer,intent(in):: tZ   ! target Z
      record /ptcl/ a(*)        ! ouptut. produced ptcl
      integer,intent(out):: ntp ! # of produced ptcls

      real(8):: TargetA, xs
      if( HowPhotoP > 0 .and.  plab.fm.p(4) > MinPhotoProdE ) then
!         call ep2cosCond
!         call cfixModel( pj )
!         call cgpXsec(Media(MediaNo).A,  E, xs)
!         prob = xs*Media(MediaNo).mbtoPX0
!         call rndc(u)
!         tgp = -log(u)/prob
!　　　−ーーーーーーーーーーー above is probably not needed

!////         call ep2cosPtcl(plab)
!         call epfixTarget2(ActiveMdl, Media(MediaNo))
!         call ep2cosCond2(Media(MediaNo).colA, 
!       *         Media(MediaNo).colZ, Media(MediaNo).colXs)
         TargetA= tA
         call cgpXsec(TargetA,  plab.fm.p(4), xs)
!////         call ep2cosCond2(tA, tZ, xs)   ! xs will not be used
         call cphotop        ! Cosmos function                     
!!         call eppushPtcl(cTrack)  ! puch Cosmos made  ptcl into Epics
         a(1:Nproduced) = Pwork(1:Nproduced)
         ntp = Nproduced
      else
         write(0,*) ' HowPhotoP=', HowPhotoP, ' Eg=', plab.fm.p(4)
         write(0,*) ' Either of above NG for photon projectile'
         stop
      endif

      end

      subroutine cutbyangle(a, ntp0,  ntp)
      implicit none
#include  "Zptcl.h"
#include  "Zcode.h"
#include  "Zevhnv.h"
#include  "Zevhnp.h"
#include  "Zmanagerp.h"
      include "Zprivate.h"
      record /ptcl/  a(*)
      integer ntp0 ! input. number of ptcls. in a
      integer ntp  ! output. could be the same as ntp0
      integer j 
      integer  i
      real*8 p, wz
      j = 0
      do i = 1, ntp0
         p = a(i).fm.p(1)**2 + a(i).fm.p(2)**2 +
     *       a(i).fm.p(3)**2
         p = sqrt(p)
         wz = a(i).fm.p(3)/p 
         if( wz .ge. wzmin .and. wz .le. wzmax ) then
            j = j + 1
            a(j)=a(i)
         endif
      enddo
      ntp = j
      end
      subroutine sortbyke(a, ntp)
      implicit none
#include  "Zptcl.h"
#include  "Zcode.h"

      include "Zprivate.h"
      integer  ntp
      record /ptcl/  a(*)
c             projectile and target information (both befor
c             and after collision ) in different system.
      integer  i
      do i = 1, ntp
         ke(i) = a(i).fm.p(4) - a(i).mass
      enddo
      call kqsortd(ke, indx, ntp)
      call ksortinv(indx, ntp)  
c       ke( indx(1) ) is the highest energy
      end
      subroutine outtrace(nev, a, ntp)
      implicit none
#include  "Zptcl.h"
#include  "Zcode.h"
#include  "Zevhnp.h"
#include  "Zevhnv.h"
#include  "Zmass.h"
#include  "Zmanagerp.h"
#include  "Zmanager.h"
#include  "Ztrackp.h"
      include  "Zprivate.h"

      integer ntp, nev
      record /ptcl/ a(ntp)
      integer  i, j, leng, icon, klena
      real*8  p
c     , wx, wy, wz 
c      real  x1, y1, z1, x2, y2, z2
      character*100  tracefile
      real(8):: x1(3), x2(3)

      write(tracefile, *) TraceDir(1:klena(TraceDir))//'/trace', nev
      call kseblk(tracefile, ' ', leng)
      call copenfw(TraceDev, tracefile(1:leng), icon)
      if(icon .ne. 0) then
         call cerrorMsg('tracefile could not be opened',0)
      endif
      if(xyz .eq. 0) then
         x1(1:3)=0.
      else
         x1(1:3)=xpos(1:3)
      endif

!          colliding partilces
      if(Trace == 1 ) then
         call puttrace(x1, pj, -2*trackl)
         call puttrace(x1, tg, -2*trackl) 
      endif

      do j = 1, ntp
         i = indx(j)
         call puttrace(x1, a(i), trackl)
      enddo
      close(TraceDev)
      end
      subroutine puttrace(x1, a, leng)
      implicit none
#include  "Zmanagerp.h"
#include  "Zptcl.h"
#include  "Ztrackp.h"
      include  "Zprivate.h"

      real(8),intent(in)::x1(3)  !   origin
      record /ptcl/ a   ! a partilce
      real(8),intent(in):: leng    ! length of track to be drawn

      real(8)::wx(3), x2(3), p

      p = sqrt( dot_product(a.fm.p(1:3), a.fm.p(1:3)))
      wx(1:3)= a.fm.p(1:3)/p
      x2(1:3) = x1(1:3) + wx(1:3)*leng
      write(TraceDev,'(3g14.5, i3, g14.4, i3, i2)') 
     *     x1(1:3),
     *     a.code,  a.fm.p(4) - a.mass, a.charge,
     *     0
      write(TraceDev, '(3g14.5, i3, g14.4, i3, g14.4)' )
     *     x2(1:3),
     *     a.code,  a.fm.p(4) - a.mass, a.charge,
     *     trackl 
      write(TraceDev, *) 
      write(TraceDev, *) 
      end

      subroutine outresulB(a, ntp)
      implicit none
#include  "Zptcl.h"
#include  "Zcode.h"
#include  "Zevhnp.h"
#include  "Zevhnv.h"
#include  "Zmass.h"
#include  "Zmanagerp.h"
#include  "Zmanager.h"
#include  "Ztrackp.h"
      include  "Zprivate.h"
c general process information; only for dpmjet3
      INTEGER IPROCE,IDNODF,IDIFR1,IDIFR2,IDDPOM,IPRON
      COMMON /POPRCS/ IPROCE,IDNODF,IDIFR1,IDIFR2,IDDPOM,IPRON(15,4)
c      IPROCE
c             1 non-diffractive inelastic
c             2 elestic 
c             3 quasi elestic vector meson prod. (photon)
c             4 central diffraction
c             5 single diff. ptcl 1
c             6 //           ptcl 2
c             7 double diff. 
c             8 direct photo-hadron
c For moore detail, see manual in Documents/CosmicRays/phojetShort.pdf
c               say, IDIFR1 classifies IPROCE=5

      integer ntp
      record /ptcl/ a(ntp)
      integer  i, j

      real*8  p, wx, wy, wz, pt, y, eta
      integer npzm, npzp, Ncht, Nchpzm, Nchpzp
      integer nw, diffcode(20) 
c//////////////
      record /ptcl/ inlab
      real(8):: ylab, etalab
      integer icon
      record /ptcl/ pcms, pl
      record /fmom/ gc
c/////////////
      npzm = 0
      npzp = 0
      Nchpzm = 0
      Nchpzp = 0

      do j = 1, ntp
         i = indx(j)
c         pt = sqrt( a(i).fm.p(1)**2 + a(i).fm.p(2)**2)
c         p= sqrt( pt**2  +    a(i).fm.p(3)**2 )
         if( a(i).fm.p(3) .gt. 0.) then
            npzp = npzp + 1
            if(a(i).charge .ne. 0) then
               Nchpzp = Nchpzp +1
            endif
         else
            npzm = npzm + 1
            if(a(i).charge .ne. 0) then
               Nchpzm = Nchpzm + 1
            endif
         endif
      enddo
      Ncht = Nchpzm + Nchpzp
      call getDiffCode(nw, diffcode)
      write(*,'("h ",i3, 6i6)' ) 
     *   diffcode(1), ntp, npzm, npzp, Ncht, Nchpzm, Nchpzp

      do j = 1, ntp
         i = indx(j)
         pt = sqrt( a(i).fm.p(1)**2 + a(i).fm.p(2)**2)
c         p= sqrt( pt**2  +    a(i).fm.p(3)**2 )
c
c         wx = a(i).fm.p(1)/p               
c         wy = a(i).fm.p(2)/p               
c         wz = a(i).fm.p(3)/p               
c         if(xyz .eq. 0) then
c            write(*,'(3i3,g14.5,3f17.13, i5)')
c     *        a(i).code, a(i).subcode, a(i).charge,
c     *        a(i).fm.p(4)-a(i).mass, wx, wy, wz, j
c         else
c            write(*,'(3i3,g14.5,1p3E11.3,0p3f17.13, i5)')
c     *        a(i).code, a(i).subcode, a(i).charge,
c     *        a(i).fm.p(4)-a(i).mass, xpos, ypos, zpos, 
c     *        wx, wy, wz, j
c         
c         endif
         call cyeta(a(i),  y, eta)
c////////////
cc           boost to target rest system
         call cbst1(1, tg, a(i), inlab)
         call cyeta(inlab,  ylab, etalab)
c///////////////
c/////////////
c         call clorez(gc, a(i), pl)
c         write(*,'(2i4, 1p, g13.4, 0p, f10.4, 1p,6g13.4)')
c     *   a(i).code,  a(i).charge,  pt,  eta, pl.fm.p(4)-pl.mass,
c     *     a(i).fm.p(1:4)
c//////////////////
         write(*,'(3i4, 1p, g13.4, 0p, f10.4,1p, 6g13.4)')
     *   a(i).code,  a(i).subcode,  a(i).charge,
     *     pt, y,  eta, a(i).fm.p(4)-a(i).mass,  a(i).fm.p(3),
     *     ylab, etalab
      enddo
      if(ntp .gt. 0 .or. outzero .eq. 1) then
         write(*, *) 
      endif
      end
