subroutine cintePhoton
      implicit none
#include  "Ztrack.h"
#include  "Ztrackv.h"
      integer i
      character*70 msg

c///////////////////
c              to see possibility of photon projectile 
c              with Primary is 'gamma' ;    see last
c              part of this file.
c      call cpredpmjet( MovedTrack.p,  Pwork, Nproduced)
c//////////

      if(IntInfArray(ProcessNo).process .eq. 'pair') then
         call cpair
      elseif(IntInfArray(ProcessNo).process .eq. 'compt') then
         call ccompt
      elseif(IntInfArray(ProcessNo).process .eq. 'photoe') then
         call cphotoEE  
      elseif(IntInfArray(ProcessNo).process .eq. 'photop') then
         call cphotop
      elseif(IntInfArray(ProcessNo).process .eq. 'cohs') then
         call ccohs
      elseif(IntInfArray(ProcessNo).process .eq. 'mpair') then
         call cmpair
      else
         write(msg, *) ' process for photon; Proc#=',
     *        ProcessNo,
     *        IntInfArray(ProcessNo).process, ' undef'
         call cerrorMsg(msg,  1)
         write(0, '("energy =",g12.4, "where=",i4, " w=",g12.3)') 
     *    MovedTrack.p.fm.p(4), MovedTrack.where, MovedTrack.wgt
         write(0, * ) " coszenith=", MovedTrack.vec.coszenith
         write(0, *) ' MoveStat=',MoveStat, 'No Of inte=',
     *       NumberOfInte
         do i = 1, NumberOfInte
            write(0,*) IntInfArray(i).process, ' dt=',
     *         IntInfArray(i).thickness,IntInfArray(i).length
         enddo
         write(0,*) 
         write(0, * ) " (dep,h)B//A==", TrackBefMove.pos.depth,
     *            TrackBefMove.pos.height,
     *            MovedTrack.pos.depth,  MovedTrack.pos.height
         stop
      endif
      end
c     *******************
      subroutine cpair
c     *******************
      implicit none

#include  "Zcode.h"
#include  "Zmass.h"
#include  "Ztrackp.h"
#include  "Ztrack.h"
#include  "Ztrackv.h"

c
      real*8 e1, e2, u, eg, cs, sn
      real*8 teta, teta1, teta2, cos1, sin1, cos2, sin2
      integer ica
      real*8 den, cvh2den
      record /track/aTrack
      record /coord/ dc, dce
      real*8 temp
c     
      eg = MovedTrack.p.fm.p(4)
      if(LpmEffect .and. eg .gt. LpmPairEmin) then
c         den = cthick2den(TrackBefMove.pos.depth)
         den = cvh2den(TrackBefMove.pos.height)  ! better
         call cpairErgLPM(eg, den, e1)
      else
         call cpairEnergy(eg, e1)
      endif

      e2=eg - e1
      if( e1 .gt. e2) then
c          store higher energy ptcl later
         temp = e1
         e1= e2
         e2 = temp
      endif
c         now e1 < e2
      
c            assign charge for e1
      call rndc(u)
      if(u .gt. .5) then
         ica=1
      else
         ica=-1
      endif
c     
      aTrack = MovedTrack
c      if(eg .lt. 100.e-3) then
c          take  pair angle if < 100 MeV
      if(eg .lt. 10) then
c          take  pair angle if < 10GeV
          call cPairAng(e1, masele, teta1) ! teta1 < pi/2
          if(teta1 .lt. 0.03d0) then
             cos1 = 1. - teta1**2/2
             sin1 = teta1
          else
             cos1 = cos(teta1)
             sin1 = sin(teta1)
          endif
c
          sin2 = sin1 * sqrt(  (e1**2-masele**2)/(e2**2-masele**2) )
          if(sin2 .lt. 0.03d0) then
             cos2 = 1.- sin2**2/2
          else
             cos2 = sqrt(1.d0 - sin2**2)
          endif
          call kcossn(cs, sn)

c         teta = masele/eg
c         teta1=teta* e2/eg
c         teta2=teta* e1/eg
c         cos1=1. - teta1**2/2
c         cos2=1. - teta2**2/2
c               sample direction cos. of 1st
c         sin1=teta1
c         sin2=teta2
c
         dc.r(1) = cs * sin1
         dc.r(2) = sn * sin1
         dc.r(3)=  cos1
         call ctransVectZ(MovedTrack.vec.w, dc, dce)
         call cmkptc(kelec, 0, ica, aTrack.p)
         aTrack.p.fm.p(4) = e1
         call csetDirCos(dce, aTrack)
         call ce2p(aTrack)
         Nproduced = Nproduced + 1
         Pwork(Nproduced) = aTrack.p
c            another electron (higher E)

         dc.r(1) = -cs*sin2
         dc.r(2) = -sn*sin2
         dc.r(3) = cos2
         call ctransVectZ(MovedTrack.vec.w, dc, dce)
         aTrack.p.fm.p(4) = e2
         call cmkptc(kelec, 0, -ica,  aTrack.p)
         call csetDirCos(dce, aTrack)
         call ce2p(aTrack)
         Nproduced = Nproduced + 1
         Pwork(Nproduced) = aTrack.p
      else
c          neglect pair angle
         aTrack.p.fm.p(4) = e1
         call ce2p(aTrack)
         call cmkptc(kelec, 0, ica, aTrack.p)
         Nproduced = Nproduced + 1
         Pwork(Nproduced) = aTrack.p
c     
         aTrack.p.fm.p(4) = e2
         call ce2p(aTrack)
         call cmkptc(kelec, 0, -ica, aTrack.p)
         Nproduced = Nproduced + 1
         Pwork(Nproduced) = aTrack.p
      endif
      end
c     ***********
      subroutine ccompt
c     ***********
      implicit none
c----      include '../Particle/Zcode.h'
#include  "Zcode.h"
c----      include 'Ztrack.h'
#include  "Ztrack.h"
c----      include 'Ztrackv.h'
#include  "Ztrackv.h"
c
      record /track/aTrack
      real*8 eg, e1, cs, sn, cosg, cose
      real*8 sine, tmp, sing
      record /coord/ dc, dce, dcg
c  
      call ccomptea(MovedTrack.p.fm.p(4), eg, e1, cosg, cose)

      call kcossn(cs,sn)
c           electron direction
      tmp=max(1.d0-cose*cose, 0.d0)
      sine=sqrt(tmp)
      dc.r(1)=cs*sine
      dc.r(2)=sn*sine
      dc.r(3)=cose
      call ctransVectZ(MovedTrack.vec.w, dc, dce)
      aTrack = MovedTrack
      call cmkptc(kelec, 0, -1, aTrack.p)
      aTrack.p.fm.p(4) = e1
      call csetDirCos(dce, aTrack)
      call ce2p(aTrack)
      Nproduced = Nproduced + 1
      Pwork(Nproduced) = aTrack.p
c            gamma dicrection
      tmp=max(1.d0-cosg*cosg, 0.d0)
      sing=sqrt(tmp)
      dc.r(1) = -cs*sing
      dc.r(2) = -sn*sing
      dc.r(3) = cosg
      call ctransVectZ(MovedTrack.vec.w, dc, dcg)
      aTrack.p.fm.p(4) = eg
      call cmkptc(kphoton, kcasg, 0, aTrack.p)
      call csetDirCos(dcg, aTrack)
      call ce2p(aTrack)
      Nproduced = Nproduced + 1
      Pwork(Nproduced) = aTrack.p


      end
c     ***********
      subroutine cphotop
c     ***********
      implicit none
#include  "Zcode.h"
#include  "Ztrack.h"
#include  "Ztrackv.h"
#include  "Ztrackp.h"
c
c      record /track/aTrack
c      integer i
c
c///////////////
c      real*8 sum
c      integer i
c      record /ptcl/work(3000)
c////////////
      integer ngen
      Nproduced = 0
c//////////////
c      write(0,*) ' in cphotop; bef cgpHad MoveE=',
c     *     MovedTrack.p.fm.p(4)
c//////////
      if( HowPhotoP == 1 ) then
         call csofia( TargetNucleonNo, TargetProtonNo, 
     *      MovedTrack.p, Pwork(Nproduced+1), ngen)
      elseif( HowPhotoP == 2 .or.
     *  ( HowPhotoP == 3 .and. MovedTrack.p.fm.p(4) < 2.5 ) ) then
         call cgpHad(TargetNucleonNo, TargetProtonNo, 
     *   MovedTrack.p, 2, Pwork(Nproduced+1), ngen)
      elseif( HowPhotoP == 4 .and. MovedTrack.p.fm.p(4) > 2.5 ) then
         call cgpHad(TargetNucleonNo, TargetProtonNo, 
     *   MovedTrack.p, 2, Pwork(Nproduced+1), ngen)
      else
         ngen = 0
      endif


c
c      call cgpHad(TargetNucleonNo, TargetProtonNo, 
c     *   MovedTrack.p, 2, work, ngen)
cc      work(1) = MovedTrack.p
cc      work(1).fm.p(4) = MovedTrack.p.fm.p(4)/2.
cc      call cmkptc(kelec, regptcl, -1, work(1))
cc      call ce2p(work(1))
cc      work(2) = work(1)
cc      call cmkptc(kelec, antip, 1, work(2))
cc      call ce2p(work(2))
cc      ngen =2 
c///////////////////
c      write(0,*) ' cgpHad nge=',ngen,' Npro=',Nproduced
c      do i = 1, ngen
c         sum = sum + Pwork(Nproduced+i).fm.p(4)
c      enddo
c      write(0,*) ' Eg=',MovedTrack.p.fm.p(4), 'out=',sum
cc      if(ngen .ne. 2 .or. Nproduced .ne. 0) then
cc         write(0, *) ' nge=',ngen,' Nproduced=',Nproduced
cc         write(0, *) ' proc=',IntInfArray(ProcessNo).process
cc         write(0,*) ' parent=',MovedTrack.p.fm.p(4),
cc     *           ' E1+E2=', work(1).fm.p(4)+work(2).fm.p(4)
cc         write(0,*) ' muon E=',Pwork(1).fm.p(4)
cc      endif
c      do i = 1, ngen
c         Pwork(i+Nproduced) = work(i)
c      enddo
c/////////////////////
      Nproduced = Nproduced + ngen
         
c     
c      aTrack = MovedTrack
c      do i = 1, Nproduced
c         aTrack.p = Pwork(i)
c         call cpush(aTrack)
c      enddo
      end
c     ****************
      subroutine ccohs
c     ****************
c        coherent scattering
      implicit none
#include  "Zcode.h"
#include  "Ztrack.h"
#include  "Ztrackv.h"
c     ******************
c        since coherent scattering is effective at
c        low energies where angular distribution can be
c        approximated by (1+cos^2) dcos, we simply use
   
      record /track/aTrack
      record /coord/ w
      real*8 cosg, tmp, cs, sn, sing
c             sample scattering angle from (1+cos^2)dcos
      call ksampRSA(cosg)
      tmp=1.d0-cosg*cosg
      sing = sqrt(tmp)
      call kcossn(cs,sn)
      
      w.r(1) = cs*sing
      w.r(2) = sn*sing
      w.r(2) = cosg
      call ctransVectZ(MovedTrack.vec.w,  w, w)
      aTrack = MovedTrack
c        energy unchaged;
      call csetDirCos(w, aTrack)
      call ce2p(aTrack)
      Nproduced = Nproduced + 1
      Pwork(Nproduced) = aTrack.p
      end
c     ****************
      subroutine cphotoEE
c     ****************
c        photo electric effect
      implicit none
#include  "Zcode.h"
#include  "Zmass.h"
#include  "Ztrack.h"
#include  "Ztrackv.h"

c     ******************
      record /track/aTrack
      record /coord/ w
      real*8  cs, sn, sing

      real*8 eout, eg, rEg, rEe, cost, a, tmp
c          essentially no shell down to 1keV.
      eg = MovedTrack.p.fm.p(4)
      eout = masele + eg
      
      aTrack = MovedTrack
      rEg=eg/masele
      rEe=(eout-masele) /masele
      if(rEe .le. 0.) then
         cost=1.
      else
         a = ( TargetAtomicN/137.0/137.0 + 2.0*rEe + rEg**2 )/
     *    (2.0*rEg*sqrt(2.0*rEe) )
         call ksampPEang(a, cost)
      endif
      
      tmp=1.d0-cost*cost
      sing = sqrt(tmp)
      call kcossn(cs,sn)
      
      w.r(1) = cs*sing
      w.r(2) = sn*sing
      w.r(2) = cost
      call ctransVectZ(MovedTrack.vec.w,  w, w)
      aTrack.p.fm.p(4) = eout
      call cmkptc(kelec, 0, -1, aTrack.p)
      call csetDirCos(w, aTrack)
      call ce2p(aTrack)
      Nproduced = Nproduced + 1
      Pwork(Nproduced) = aTrack.p
      end
c     *******************
      subroutine cmpair
c         magnetic pair creation
c     *******************
      implicit none
#include  "Zcode.h"
#include  "Ztrack.h"
#include  "Ztrackv.h"
c
      real*8 e1, e2, u
      integer ica, nc
      record /track/aTrack
c       
      call cmPairE(Xai, e2, nc)
c           e2 is higher energy fraction; change to real energy
c         store later in  working array, then higher one is
c        stored first in the stack to save the memory.
      e2 = MovedTrack.p.fm.p(4) * e2
      e1=MovedTrack.p.fm.p(4) - e2
c            assign charge for e1
      call rndc(u)
      if(u .gt. .5) then
         ica=1
      else
         ica=-1
      endif
c     
      aTrack = MovedTrack
      aTrack.p.fm.p(4) = e1
      call ce2p(aTrack)
      call cmkptc(kelec, 0, ica, aTrack.p)
      Nproduced = Nproduced + 1
      Pwork(Nproduced) = aTrack.p

      aTrack.p.fm.p(4) = e2
      call ce2p(aTrack)
      call cmkptc(kelec, 0, -ica, aTrack.p)
      Nproduced = Nproduced + 1
      Pwork(Nproduced) = aTrack.p
      end
      
c/////////////
      subroutine cpredpmjet(pj,  a, ntp)
#include  "Zair.h"
#include "Zptcl.h"
      record /ptcl/  pj
      integer ntp
      record /ptcl/ a(*)
      TargetNucleonNo=14
      TargetProtonNo= 7 
      call cdpmjet( pj, TargetNucleonNo, TargetProtonNo,
     *   a, ntp)
      write(0,*) ' ntp =',ntp
      stop
      end
c/////////////