1 jones 1.1       subroutine guhadr
  2                 implicit none
  3           c
  4           c*keep,gctrak.
  5           *-- author :
  6                 common/gctrak/vect(7),getot,gekin,vout(7),nmec,lmec(30),namec(30)
  7                + ,nstep ,maxnst,destep,destel,safety,sleng ,step  ,snext ,sfield
  8                + ,tofg  ,gekrat,upwght,ignext,inwvol,istop ,idecad,iekbin
  9                + , ilosl, imull,ingoto,nldown,nlevin,nlvsav,istory
 10                 real*8 pout,px,py,pz
 11                 real*8 xhatf,yhatf,zhatf,xhatr,yhatr,zhatr
 12                 real*8 xhat1f,yhat1f,zhat1f,xhat1r,yhat1r,zhat1r
 13                 real*8 psif,tthetaf,thetasc,phisc,psisc
 14                 real*8 theta_sc,phi_sc,psi_sc,x_sc,y_sc,z_sc,r_sc
 15                 real*8 phi_az,theta_az,pi,anorm
 16                 real*8 a,b,c,d,ac,p,r,pp
 17                 real*8 term1,term2,term3,term4
 18                 real*8 prob,pmax,ratio,pmax1,pmax2
 19           c
 20                 real  vect,getot,gekin,vout,destep,destel,safety,sleng ,step
 21                + ,snext,sfield,tofg  ,gekrat,upwght,iii,ii
 22 jones 1.1 c end gctrak
 23                 real rndm(3)
 24                 integer nmec,lmec,namec,nstep ,maxnst,ignext,inwvol,istop
 25                + ,idecad,iekbin,ilosl, imull,ingoto,nldown,nlevin,nlvsav,istory
 26           
 27           c      include 'geant_local.h'
 28           
 29                 pi=3.14159265
 30           
 31           c      write(*,*)'Calling hadronic interaction package ...'
 32           c      pout=vect(7)*1000.
 33           c      px=pout*vect(4)
 34           c      py=pout*vect(5)
 35           c      pz=pout*vect(6)
 36           c      xhatf=px/pout
 37           c      yhatf=py/pout
 38           c      zhatf=pz/pout
 39           c      anorm=sqrt(xhatf**2+yhatf**2+zhatf**2)
 40           c      xhatf=xhatf/anorm
 41           c      yhatf=yhatf/anorm
 42           c      zhatf=zhatf/anorm
 43 jones 1.1 c     write(*,*)xhatf,yhatf,zhatf,sqrt(xhatf**2+yhatf**2+zhatf**2)
 44           
 45           c      psif=dacos(sqrt(yhatf**2+zhatf**2))
 46           c      if(xhatf.lt.0.0) psif=-1.0*psif
 47           c      tthetaf=yhatf/zhatf
 48           
 49           c      xhat1f=xhatf*dcos(psif)-
 50           c     &     yhatf*dsin(datan(tthetaf))*dsin(psif)-
 51           c     &     zhatf*dcos(datan(tthetaf))*dsin(psif)
 52           c      yhat1f=
 53           c     &     yhatf*dcos(datan(tthetaf))-
 54           c     &     zhatf*dsin(datan(tthetaf))
 55           c      zhat1f=xhatf*dsin(psif)+
 56           c     &     yhatf*dsin(datan(tthetaf))*dcos(psif)+
 57           c     &     zhatf*dcos(datan(tthetaf))*dcos(psif)
 58           c      write(*,*)xhat1f,yhat1f,zhat1f,sqrt(xhat1f**2+yhat1f**2+zhat1f**2)
 59           
 60           c      pp=pout/1000.0-1.4
 61           c      r=pout/1000.0
 62                 
 63            100  format(1x,3(f8.3,1x))
 64 jones 1.1 c
 65           8876  call gheish
 66           c
 67           c      pout=vect(7)*1000.
 68           c      px=pout*vect(4)
 69           c      py=pout*vect(5)
 70           c      pz=pout*vect(6)
 71           c      xhatr=px/pout
 72           c      yhatr=py/pout
 73           c      zhatr=pz/pout
 74           c      anorm=sqrt(xhatr**2+yhatr**2+zhatr**2)
 75           c      xhatr=xhatr/anorm
 76           c      yhatr=yhatr/anorm
 77           c      zhatr=zhatr/anorm
 78           
 79           c      xhat1r=xhatr*dcos(psif)-
 80           c     &     yhatr*dsin(datan(tthetaf))*dsin(psif)-
 81           c     &     zhatr*dcos(datan(tthetaf))*dsin(psif)
 82           c      yhat1r=
 83           c     &     yhatr*dcos(datan(tthetaf))-
 84           c     &     zhatr*dsin(datan(tthetaf))
 85 jones 1.1 c      zhat1r=xhatr*dsin(psif)+
 86           c     &     yhatr*dsin(datan(tthetaf))*dcos(psif)+
 87           c     &     zhatr*dcos(datan(tthetaf))*dcos(psif)
 88           c      write(*,*)xhat1r,yhat1r,zhat1r
 89           c
 90            
 91           c      psisc=dasin(xhat1r)
 92           c      if(dcos(psisc).ne.0.0D0) then 
 93           c         thetasc=dasin(yhat1r/dcos(psisc))
 94           c      else
 95           cc         write(*,*)'Zero problem with psisc'
 96           cc         theta_az=100.0
 97           cc         phi_az=400.0
 98           c         goto 3344
 99           c      endif
100           c      if(dcos(thetasc).ne.0.0D0) then 
101           c         phisc=datan(dtan(psisc)/dcos(thetasc))
102           c      else
103           cc         write(*,*)'Zero problem with thetasc'
104           cc         theta_az=100.0
105           cc         phi_az=400.0
106 jones 1.1 cc         goto 3344
107           c      endif
108           
109           c      theta_sc=thetasc
110           c      phi_sc=phisc
111           
112           cc
113           cc Another method of changing to spherical coordinates
114           cc
115           c      psi_sc=datan(dtan(phi_sc)*dcos(theta_sc))
116           c      x_sc=dsin(psi_sc)
117           c      y_sc=dcos(psi_sc)*dsin(theta_sc)
118           c      z_sc=dcos(psi_sc)*dcos(theta_sc)
119           c      r_sc=sqrt(x_sc**2+y_sc**2)
120           c      if(z_sc.ne.0.0.and.y_sc.ne.0.0) then
121           c         theta_az=datan(r_sc/z_sc)
122           c         phi_az=datan(x_sc/y_sc)
123           c      else
124           cc         write(*,*)'Zero problem with zsc or ysc'
125           cc         theta_az=100.0
126           cc         phi_az=400.0
127 jones 1.1 c         goto 3344         
128           c      endif
129           c      theta_az=theta_az*180.0/pi
130           c      phi_az=phi_az*180.0/pi
131           c      if(y_sc.lt.0.0) phi_az=phi_az+180.0
132           c      if(y_sc.ge.0.0.and.x_sc.lt.0.0) phi_az=phi_az+360.0
133           c      geant_theta=theta_az
134           c      geant_phi=phi_az
135           c
136           c      a=acoeff(1,1)+acoeff(1,2)*pp+acoeff(1,3)*pp**2+
137           c     >   acoeff(1,4)*pp**3+acoeff(1,5)*pp**4+acoeff(1,6)*pp**5
138           c      b=acoeff(2,1)+acoeff(2,2)*pp+acoeff(2,3)*pp**2+
139           c     >   acoeff(2,4)*pp**3+acoeff(2,5)*pp**4+acoeff(2,6)*pp**5
140           c      c=acoeff(3,1)+acoeff(3,2)*pp+acoeff(3,3)*pp**2+
141           c     >   acoeff(3,4)*pp**3+acoeff(3,5)*pp**4+acoeff(3,6)*pp**5
142           c      d=acoeff(4,1)+acoeff(4,2)*pp+acoeff(4,3)*pp**2+
143           c     >   acoeff(4,4)*pp**3+acoeff(4,5)*pp**4+acoeff(4,6)*pp**5
144           
145           c      p=r
146           c      r=r*dsin(theta_az/180.0*pi)
147           c      term1=a*r
148 jones 1.1 c      term2=b*r**2
149           c      term3=c*r**4
150           c      term4=d*p*dsin(5.0*theta_az/180.0*pi)
151           c      ac=term1/(1+term2+term3)+term4
152           c
153           cc Now we have the values of phi, ac(theta) and pn,pt.  Now calculate
154           cc the relative probability.
155           cc
156           c      pmax1=1+pn*ac*dcos(datan(+1.0*pt/pn))+
157           c     >     pt*ac*dsin(datan(+1.0*pt/pn))
158           c      pmax2=1+pn*ac*dcos(datan(+1.0*pt/pn)+pi)+
159           c     >     pt*ac*dsin(datan(+1.0*pt/pn)+pi)
160           c      pmax=max(pmax1,pmax2)
161           c      prob=1+pn*ac*dcos(phi_az*pi/180.0)+
162           c     >	         pt*ac*dsin(phi_az*pi/180.0)
163           c      ratio=prob/pmax
164           cc      write(*,*)pn/0.384/(-0.87524),pt/(0.384)
165           c      call grndm(rndm,3)
166           cc      write(*,*)rndm(1)
167           c      if(rndm(1).gt.ratio) then
168           cc	  write(*,*)'Regeneration ...'
169 jones 1.1 c          goto 8876
170           c      endif
171           
172           c 3344 continue
173           cc      write(*,*)'Returning from guhadr ...' 
174                 return
175                 end