Gencol.f

Go to the documentation of this file.

#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
      type(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
!                 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(outresul == 0 ) then
            call outresula(w, ntp)
         elseif(outresul == 1 ) then
            call outresulb(w, ntp)
         elseif(outresul == 2 ) then
            call outresulc(w, ntp)
         elseif(outresul == 3 ) then
            call outresulx(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, outresul)
!       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)
!          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
!            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
!         call csibyllinit
!         ActiveMdl = 'sibyll'
      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
!c         call epReadMediaForMuon  !////////////
      endif
      end
      subroutine  epreadmediaformuon
#include  "Zcode.h"
#include  "Zptcl.h"
      
      include  "Zprivate.h"
      integer icon

      if( targetname /= " " ) then
         mediano = 1
!ccc         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
!         read(inpfileno,*, end=100) xpos, ypos, zpos
         read(inpfileno,*, end=100) xpos(1:3)
      endif
      eof = 0
      return
 100  continue
      eof = 1
      end
!     *******************
      subroutine formpjtg(confirm)
!     ******************
      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"
      type(ptcl)::resttg, zcmsp, zplab
      integer confirm  !  input. if 0, root s is  not printed.
                      !         else  printed
      real*8   roots, s
!         form projectile  and target
      integer::icon

      call cmkptc(pjcode, pjsub, pjchg,   pj)
!       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)
         
!       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)
!       
      if(tgpx .eq. 0. .and. tgpy .eq. 0. .and.
     *        tgpz .eq. 0.)  then
!     target is at rest;  s = (Ep+Et)^2 - (Pp+Pt)^2
!                           = (Ep+Mt)^2 - Pp^2
!                           =  Mp^2 +2EpMt +Mt^2
!   
         s= 2*pj%fm%p(4)*tg%mass +tg%mass**2 + pj%mass**2
      else
!         by  general formula;
!               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)
!      if(confirm .ne. 0) then
         write(0, *) ' roots/2=', roots/2
         write(0, *) 's,roots above are total not by /n'
!      endif
!c           boost to target rest system
      call cbst1(1, tg, pj, plab)
!         make rest target
      call cbst1(2, tg, tg, resttg)
       !   get /n cms
      call cgeqm2(plab, resttg, cmsp, icon)   
!////////////
      write(0,*) ' Cmsp= (/n)',cmsp%fm%p(:), cmsp%mass
      write(0,*) ' plab=',plab%fm%p(4)
      write(0,*) ' tg=', tg%fm%p(4)
!///////////
      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)   
!         inform Zcmsp to cintModels.f 
!      call cputGencolCMS(Cmsp) 
      call cputgencolcms(zcmsp) 
      end
!     ************************
      subroutine outresula(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
      type(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)
         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
         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


         else
            write(*,'(i3,i5,i4,g14.5,1p3E11.3,0p3f17.13, i8)')
     *        a(i)%code, a(i)%subcode, a(i)%charge,
     *        a(i)%fm%p(4)-a(i)%mass, xpos(1:3), wx(1:3), 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 "Zmass.h"
#include  "Zevhnv.h"
#include  "Zevhnp.h"
#include  "Zmanagerp.h"
      include "Zprivate.h"

      type(ptcl)::  a(*)
!             projectile and target information (both befor
!             and after collision ) in different system.
!
      integer  ntp
      integer j
      integer tZ, tA, L
      real*8  xs, SIGbdif, xsair
      real*8   roots
      
      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
!!!!!!!         
!     you may include  ToPrintXS.f here to print out E vs XS for Sibyll
 !!!!!!!!
!           generate event
         if(activemdl .eq. 'qgsjet1') then
#ifdef QGSJET1
            call qgs01event(plab, ta, tz, a, ntp)
#endif
!!         elseif(ActiveMdl .eq. 'sibyll') then
!!!
!!            call csibyllevent(plab, tA, tZ, a, ntp)
         else
            call chacol(plab, ta, tz,  a, ntp)
         endif
      endif
      call getimpactparam(b)
      do j = 1, ntp
!               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
!  #include  "Zcode.h"
#include  "Ztrack.h"
#include  "Ztrackv.h"
#include  "Ztrackp.h"
#include  "Zelemagp.h"

!       include "Zprivate.h"

      type(ptcl):: plab  ! input proj. in Lab
      integer,intent(in):: tA   ! target A
      integer,intent(in):: tZ   ! target Z
      type(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"
      type(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
      type(ptcl)::  a(*)
!             projectile and target information (both befor
!             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)  
!       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
      type(ptcl):: a(ntp)
      integer  i, j, leng, icon, klena
      real*8  p
!     , wx, wy, wz 
!      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
      type(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)))
      if( p >0. ) then
         wx(1:3)= a%fm%p(1:3)/p
      else
         wx(1:3) = 0.
      endif
      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"
      integer,save:: nthev=0
      include  "Zprivate.h"
! general process information; only for dpmjet3
      INTEGER IPROCE,IDNODF,IDIFR1,IDIFR2,IDDPOM,IPRON
      COMMON /poprcs/ iproce,idnodf,idifr1,idifr2,iddpom,ipron(15,4)
!      IPROCE
!             1 non-diffractive inelastic
!             2 elestic 
!             3 quasi elestic vector meson prod. (photon)
!             4 central diffraction
!             5 single diff. ptcl 1
!             6 //           ptcl 2
!             7 double diff. 
!             8 direct photo-hadron
! For moore detail, see manual in Documents/CosmicRays/phojetShort.pdf
!               say, IDIFR1 classifies IPROCE=5

      integer ntp
      type(ptcl):: a(ntp)
      integer  i, j

      integer:: npp
      

      
      real*8  p, wx, wy, wz, pt, y, eta
      integer npzm, npzp, Ncht, Nchpzm, Nchpzp
      integer nw, diffcode(20)
      real(8):: Xfpz, Xfke
!//////////////
      type(ptcl):: inlab
      real(8):: ylab, etalab
      integer icon
      type(ptcl):: pcms, pl
      type(fmom):: gc
!/////////////
      nthev=nthev+1
      
      npzm = 0
      npzp = 0
      nchpzm = 0
      nchpzp = 0

      
      do j = 1, ntp
         i = indx(j)
         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
!      ntpp = 0   !  ?

      do j = 1, ntp
         i = indx(j)
         pt = sqrt( a(i)%fm%p(1)**2 + a(i)%fm%p(2)**2)
         p= sqrt( pt**2  +    a(i)%fm%p(3)**2 )
         xfke = ( a(i)%fm%p(4)-a(i)%mass) /( pj%fm%p(4)- pj%mass)
         xfpz =  a(i)%fm%p(3) / pj%fm%p(3)
         
         wz = a(i)%fm%p(3)/p               
         call cyeta(a(i),  y, eta)
         write(*,'(3i4, 1p, 7g13.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),
     *        xfke, xfpz           
!     *     ylab, etalab
!         write(*,'(3i4, 1p, 3g13.4)')
!     *   a(i)%code,  a(i)%subcode,  a(i)%charge,
!     *    a(i)%fm%p(4)-a(i)%mass,    a(i)%fm%p(3), wz

         
      enddo
      if(ntp .gt. 0 .or. outzero .eq. 1) then
         write(*, *) 
      endif
      end
      subroutine outresulc(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"
      integer,save:: nthev=0
      include  "Zprivate.h"
! general process information; only for dpmjet3
      INTEGER IPROCE,IDNODF,IDIFR1,IDIFR2,IDDPOM,IPRON
      COMMON /poprcs/ iproce,idnodf,idifr1,idifr2,iddpom,ipron(15,4)
!      IPROCE
!             1 non-diffractive inelastic
!             2 elestic 
!             3 quasi elestic vector meson prod. (photon)
!             4 central diffraction
!             5 single diff. ptcl 1
!             6 //           ptcl 2
!             7 double diff. 
!             8 direct photo-hadron
! For moore detail, see manual in Documents/CosmicRays/phojetShort.pdf
!               say, IDIFR1 classifies IPROCE=5

      integer ntp
      type(ptcl):: a(ntp)
      integer  i, j, icon

      
      integer:: npp
      
      real*8  p, wx, wy, wz, pt, y, eta
      integer npzm, npzp, Ncht, Nchpzm, Nchpzp
      integer nw, diffcode(20)
      real(8):: Xfpz, Xfke
!//////////////
      type(ptcl):: restn
      type(ptcl):: cmptcl
      type(fmom):: gc
      real(8):: s, rootsh, xf, xflab, xfmax
      real(8):: Esmax, Pszmax, Pzmax, Ebypzmax, Esbyrootsh
      real(8):: ys, etas, Es, Elab, pabs, cost, teta
!/////////////
      nthev=nthev+1


      call cmkptc(6, -1, 1, restn)
      restn%fm%p(1:3)=(/0.,0.,0./)
      restn%fm%p(4) =  restn%mass
      call cgetcm(pj, restn, gc, icon)
      s =( restn%mass + pj%fm%p(4) ) *2 * restn%mass
      rootsh = sqrt(s)/2
      esmax=rootsh - pj%mass
      pszmax=sqrt( esmax**2 - maspic**2 )    ! pion case
      pzmax = gc%p(4)*pszmax + gc%p(3)* esmax !

      do j = 1, ntp
         call cbst0(j, gc, a(j), cmptcl) ! j==1  must be called.
         call cyeta(cmptcl,  ys, etas)
         call cyeta(a(j), y, eta)

         elab = a(j)%fm%p(4)
         es = cmptcl%fm%p(4)
         xf = cmptcl%fm%p(3)/rootsh
         xflab = a(j)%fm%p(3)/pzmax
         pt =sqrt( a(j)%fm%p(1)**2 + a(j)%fm%p(2)**2 )
         pabs = sqrt( pt**2 + a(j)%fm%p(3)**2 )
         cost = a(j)%fm%p(3)/pabs
         teta = acos(cost)
         write(*, '(2i4, 1p, 9g14.4)')
     *     a(j)%code,  a(j)%charge, etas, es, xf, eta, elab,
     *     xflab,  cost, teta, pt

      enddo
      if(ntp .gt. 0 .or. outzero .eq. 1) then
         write(*, *) 
      endif
      end


      subroutine outresulx(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
      type(ptcl):: a(ntp)
      integer  i, j
      real*8  p, mass2, p2
      real*8  wx(3)
      integer nw, difcode(20)

!      call getDiffCode(nw, dffcode)
!!!
      call csibylgetdiffcode(nw, difcode)
!      write(*,'("h ",i3, 7i6)' )
!     * difcode(1),ntp, npzm, npzp, Ncht, Nchpzm, Nchpzp, nthev
!      
!      do i = 1, nw
!         write(*,'(a, 3i10)')
!     *         'dif ', i, difcode(i), ntp
!      enddo
      write(*,*) 'dif ', nw, difcode(1:nw)
!!!!!!      
      do i = 1, ntp
         p2 =  dot_product(a(i)%fm%p(1:3), a(i)%fm%p(1:3) )
         p = sqrt(p2)
!     mass2= a(i)%fm%p(4)**2 - p2
         mass2= (a(i)%fm%p(4)- p)*(a(i)%fm%p(4)+ p)
!         if( a(i)%code  /= 1 .and.  
!     *        (  abs(a(i)%mass/sqrt(mass2) -1.) > 0.01  )
!     *           .or.
!     *        a(i)%code == 0 .and.
!     *        (  abs( a(i)%mass-sqrt(mass2) ) > 0.03  ) ) then       
!!!            write(*,
!!!     *     '(i3,i5,i4,1p, 6g22.14)')
!!!     *        a(i)%code, a(i)%subcode, a(i)%charge,
!!!     *           a(i)%fm%p(:), a(i)%mass, sqrt(mass2)
!         endif
      enddo
      if(ntp .gt. 0 .or. outzero .eq. 1) then
         write(*, *) 
      endif
      end