EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
hijels.f
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file hijels.f
1 C
2 C
3 C*******************************************************************
4 CThis subroutine performs elastic scattering between two nucleons
5 C
6 C*******************************************************************
7  SUBROUTINE hijels(PSC1,PSC2)
8  IMPLICIT DOUBLE PRECISION(d)
9  dimension psc1(5),psc2(5)
10  common/hiparnt/hipr1(100),ihpr2(50),hint1(100),ihnt2(50)
11  SAVE /hiparnt/
12  common/ranseed/nseed
13  SAVE /ranseed/
14 C
15  cc=1.0-hint1(12)/hint1(13)
16  rr=(1.0-cc)*hint1(13)/hint1(12)/(1.0-hipr1(33))-1.0
17  bb=0.5*(3.0+rr+sqrt(9.0+10.0*rr+rr**2))
18  ep=sqrt((psc1(1)-psc2(1))**2+(psc1(2)-psc2(2))**2
19  & +(psc1(3)-psc2(3))**2)
20  IF(ep.LE.0.1) RETURN
21  els0=98.0/ep+52.0*(1.0+rr)**2
22  pcm1=psc1(1)+psc2(1)
23  pcm2=psc1(2)+psc2(2)
24  pcm3=psc1(3)+psc2(3)
25  ecm=psc1(4)+psc2(4)
26  am1=psc1(5)**2
27  am2=psc2(5)**2
28  amm=ecm**2-pcm1**2-pcm2**2-pcm3**2
29  IF(amm.LE.psc1(5)+psc2(5)) RETURN
30 C ********elastic scattering only when approaching
31 C to each other
32  pmax=(amm**2+am1**2+am2**2-2.0*amm*am1-2.0*amm*am2
33  & -2.0*am1*am2)/4.0/amm
34  pmax=abs(pmax)
35 20 tt=atl_ran(nseed)*min(pmax,1.5)
36  els=98.0*exp(-2.8*tt)/ep
37  & +52.0*exp(-9.2*tt)*(1.0+rr*exp(-4.6*(bb-1.0)*tt))**2
38  IF(atl_ran(nseed).GT.els/els0) go to 20
39  phi=2.0*hipr1(40)*atl_ran(nseed)
40 C
41  dbx=pcm1/ecm
42  dby=pcm2/ecm
43  dbz=pcm3/ecm
44  db=sqrt(dbx**2+dby**2+dbz**2)
45  IF(db.GT.0.99999999d0) THEN
46  dbx=dbx*(0.99999999d0/db)
47  dby=dby*(0.99999999d0/db)
48  dbz=dbz*(0.99999999d0/db)
49  db=0.99999999d0
50  WRITE(6,*) ' (HIJELS) boost vector too large'
51 C ********Rescale boost vector if too close to unity.
52  ENDIF
53  dga=1d0/sqrt(1d0-db**2)
54 C
55  dp1=sqrt(tt)*sin(phi)
56  dp2=sqrt(tt)*cos(phi)
57  dp3=sqrt(pmax-tt)
58  dp4=sqrt(pmax+am1)
59  dbp=dbx*dp1+dby*dp2+dbz*dp3
60  dgabp=dga*(dga*dbp/(1d0+dga)+dp4)
61  psc1(1)=dp1+dgabp*dbx
62  psc1(2)=dp2+dgabp*dby
63  psc1(3)=dp3+dgabp*dbz
64  psc1(4)=dga*(dp4+dbp)
65 C
66  dp1=-sqrt(tt)*sin(phi)
67  dp2=-sqrt(tt)*cos(phi)
68  dp3=-sqrt(pmax-tt)
69  dp4=sqrt(pmax+am2)
70  dbp=dbx*dp1+dby*dp2+dbz*dp3
71  dgabp=dga*(dga*dbp/(1d0+dga)+dp4)
72  psc2(1)=dp1+dgabp*dbx
73  psc2(2)=dp2+dgabp*dby
74  psc2(3)=dp3+dgabp*dbz
75  psc2(4)=dga*(dp4+dbp)
76  RETURN
77  END