c Read file created with MARS and write for003.dat for ICOOL (many)
      implicit double precision (a-h,o-z)
c       real*8 xp,pp,pol,tp,evtwt,tbref
      character*80 probtitle
      dimension xp(3),pp(3),pol(3),xbref(3),pbref(3),xmass(12)
      data (xmass(i), i=1,12) /0.9383d0,0.93956d0,0.1396d0,0.1396d0,
     !  0.4937d0,0.4937d0,0.1057d0,0.1057d0,0.d0,0.511d-3,0.511d-3,
     !  0.9383d0/
        open(unit=3,file='E=10.GeV.txt',status='unknown',iostat=ioc)
        open(unit=4,file='pion-pos.dat',status='unknown',iostat=ioc)
        open(unit=44,file='pion-neg.dat',status='unknown',iostat=ioc)
        open(unit=20,file='proton.dat',status='unknown',iostat=ioc)
        if( ioc .ne. 0 ) go to 800
        ipstart=1
        ipstop=1000001
        tbref=0.0d0
        do i=1,3
        xbref(i)=0.0d0
        pbref(i)=0.0d0
        pol(i)=0.0d0
        end do
        pol(3)=1.d0    ! maximum polarization
        probtitle='MARS run: E=10 GeV; tgt=tantalum ON=rod'
        nfile=4
        nfileq=44
        write(nfile,'(a80)') probtitle     ! skip title card on beam input file
        write(nfile,*) xbref,pbref,tbref   ! ref particle info
        write(nfileq,'(a80)') probtitle    ! skip title card on beam input file
        write(nfileq,*) xbref,pbref,tbref   ! ref particle info
c
        do 200 ip=ipstart,ipstop    ! loop over particles
!   
        iflag = 0
        ipnum=1
        ipflg=0
!        xke is kinetic energy in GeV;
          read(3,*,iostat=ioc)ievt,iptyp,xke,evtwt,xp,pp,tp
! 
        if(ioc .ne. 0 ) then   ! eof reading FOR003 beam input
          go to 300
        end if
! transformation here
c         energy=xke+xm_pi   ! total energy 
         do i=1,3
         xp(i)=xp(i)*0.01d0    !conversion cm ===> m
c pc=sqrt(xke*xke+2.*xke*mc^2)
         factor=sqrt(xke*(xke+2.d0*xmass(iptyp)))
         pp(i)=pp(i)*factor   !!! pp is canonical momentum NOW
         end do
!         print*,'Npart',ip
! select protons and anti-protons
        if(iptyp.eq.1.or.iptyp.eq.12) then
          write(20,*,iostat=ioc)ievt,iptyp,xke,evtwt,xp,pp,tp
        go to 200
        endif
c        if(ip.le.10000)then
        if(iptyp.eq.3.or.iptyp.eq.5.or.iptyp.eq.7) then
        write(nfile,*,iostat=ioc)ievt,ipnum,iptyp,ipflg,tp,evtwt,xp,pp,
     &     pol
         else if(iptyp.eq.4.or.iptyp.eq.6.or.iptyp.eq.8) then
c          iptyp=4
        write(nfileq,*,iostat=ioc)ievt,ipnum,iptyp,ipflg,tp,evtwt,xp,pp,
     &     pol
         endif
 200     continue
        goto 900
800   continue
c  Error in input
      write(*,*)' <<<<<     ERROR OPENING FILE     >>>>>'
      write(2,*)' <<<<<     ERROR OPENING FILE     >>>>>'
      iflag = -2
      go to 900
 300     continue
      write(*,*)' <<<<<< ERROR EOF   >>>>>'
      write(2,*)' <<<<<< ERROR EOF   >>>>>'
900   continue
      stop
      end