1 jones 1.1       subroutine gukine
  2           c
  3           c  jpsullivan may 20, 1993: a simple example of gukine
  4           c
  5                 implicit none
  6                 integer       ikine,itra,istak,ivert,ipart,itrtyp,napart,ipaold
  7                 real          pkine,amass,charge,tlife,vert,pvert
  8                 common/gckine/ikine,pkine(10),itra,istak,ivert,ipart,itrtyp
  9                +      ,napart(5),amass,charge,tlife,vert(3),pvert(4),ipaold
 10                 integer nubuf
 11                 parameter (nubuf=10)
 12                 real ubuf(nubuf),rndm(3)
 13                 integer ntbeam,nttarg,nvtx
 14                 save nvtx
 15                 real plab(3),ptot,etot
 16                 real*8 ang,angr,phi,phir,psir
 17                 real mp,ea,pa
 18                 real r2,r3
 19                 integer nt,ierri,i0,i1,i2,i3
 20                 integer i,j,k,ichoice,ichoice2
 21                 real*8 rotmat,xyz(3),xyznew(3),termang
 22 jones 1.1       integer ic,ii,jj
 23                 integer*4 junk1
 24                 real*8 xfp,yfp,tthfp,tphfp,pfp,junk2,junk3,thfp,phfp
 25                 real*8 e0(10),eang(10),hang(10),targ_thick(10)
 26                 real*8 escat,pscat
 27                 integer*4 ikinsetting
 28                 real*8 fg
 29                 real dptgt
 30                 character*80 junkline
 31           c
 32                 include 'fpp_local.h'
 33                 include 'geant_local.h'
 34           c
 35           c      include 'parameter.h'
 36           c      include 'espace_type.h'
 37           c      include 'detector.h'
 38           c      include 'transport.h'
 39           c      include 'option.h'
 40           c
 41           c      common/kincom/rotmat(3,3)
 42           c
 43 jones 1.1 111   format(a80)
 44                 write(*,*)'nevent =',nevent
 45                 mp=938.2785
 46                 fg=3.14159265/180.0
 47           
 48                 if(nevent.eq.0) then
 49           c         write(*,*)'Which kinematics setting ?'
 50           	 open(unit=1,file='geant_kinematics.dat',type='UNKNOWN')
 51                    read(1,*)ikinsetting
 52                    write(*,*)'Kinematics Setting: ',ikinsetting
 53           	 close(unit=1)
 54                    open(unit=1,file='hdr_gep.dat',status='old')
 55                    do i=1,10
 56                       read(1,*)e0(i),eang(i),hang(i),targ_thick(i)
 57                       write(*,*)e0(i),eang(i),hang(i),targ_thick(i)
 58                    enddo
 59                    einc=e0(ikinsetting)
 60                    hrse_ang=eang(ikinsetting)
 61                    hrsh_ang=-1.0*hang(ikinsetting)
 62                    escat=mp/(1.0+mp/einc-cos(fg*hrse_ang))
 63                    pscat=sqrt(einc**2+escat**2-
 64 jones 1.1      $     2.0*einc*escat*cos(fg*hrse_ang))
 65                    pcentral=pscat
 66                    close(unit=1)
 67                 endif
 68           
 69                 nevent=nevent+1
 70           
 71                 call grndm ( rndm , 3 )
 72                 
 73                 dptgt = -0.05+rndm(1)*0.100
 74           c
 75           ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
 76           c
 77           c  make the vertex at z=0, (x,y) chosen randomly in a 20 cm x 20 cm square
 78           c
 79           c        write ( 6,* ) ' gukine called',icount
 80           c
 81           	vert(1)=-65.0+rndm(2)*110.0
 82           	vert(2)=-6.0+rndm(3)*12.0
 83                   vert(3) = 0.0
 84                   pfp=dptgt     
 85 jones 1.1 
 86           c        write(*,*)xfp,phfp,yfp,thfp
 87           
 88            1135   dpmom=pfp
 89           	pfp=pcentral*(1.0+pfp)
 90                   pmom=pfp
 91                   ea=sqrt(pfp**2+mp**2)
 92                   tinit=ea-mp
 93                   ntbeam  = 0.0
 94                   nttarg  = 0.0
 95                   ubuf(1) = 0.0
 96           c        write(6,*)'***********************************************'
 97           c      write ( 6,* ) ' guout:  xyz =',vert(1),vert(2),vert(3)
 98           c
 99           cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
100           c
101           1000  continue
102           c
103                 ipart      = 14       ! geant pid  (8=pi+,9,pi-, 14 =p)
104           	ptot=pfp/1000.
105                   call grndm ( rndm , 3 )
106 jones 1.1         ang=-2.2+rndm(1)*4.4
107                   phi=-7.0+rndm(2)*12.5
108                   angr=ang*3.14159265/180.0
109                   phir=phi*3.14159265/180.0
110                   psir=datan(dtan(phir)*dcos(angr))
111           c        ang=(-3.66+rndm(1)*7.32)/100.
112           c        phi=(-9.98+rndm(2)*19.96)/100.
113           c        ang=5.0
114           c        phi=5.0
115           c
116           c these are the actual parameters for the track.  
117           c
118           
119                 call gsvert ( vert,ntbeam,nttarg,ubuf,0,nvtx )
120                 etot       = ea - mp
121                 plab(1)    = ptot*dsin(psir)
122                 plab(2)    = ptot*dsin(angr)*dcos(psir)
123                 plab(3)    = ptot*dcos(angr)*dcos(psir)
124                 
125                  write(6,*)'initial momentum = ',ptot
126                  write(6,*)'initial components are',plab(1),plab(2),plab(3)
127 jones 1.1 c
128           c      write ( 6,902 ) ipart
129           c902   format ( ' **********************************',
130           c     x         ' gukine: start new track: ipart=',i5 )
131           c
132           c        write ( 6,* ) ' filling etot histo'
133                 call hfill ( 100, etot, 0., 1. )
134           c
135           c  give the track to geant
136           c
137           c        write ( 6,* ) ' gskine called'
138                 call gskine ( plab,ipart,nvtx,ubuf,0,nt )
139           c
140                 if ( nt.le.0 ) then
141                     write ( 6,* ) ' gukine: error defining track'
142                     write ( 6,* ) '         i=',i,' nt=',nt
143                     stop
144                 end if
145           c
146           c  print the vertex and track info
147           c
148 jones 1.1 c     call gpvert(0)
149           c     call gpkine(0)
150           c
151                 return
152           c
153                 end