Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
hijcsc.f
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file hijcsc.f
1 C
2 C*******************************************************************
3 C This subroutine performs elastic scatterings and possible
4 C elastic cascading within their own nuclei
5 c*******************************************************************
6  SUBROUTINE hijcsc(JP,JT)
7  dimension psc1(5),psc2(5)
8  common/hijcrdn/yp(3,300),yt(3,300)
9  SAVE /hijcrdn/
10  common/hiparnt/hipr1(100),ihpr2(50),hint1(100),ihnt2(50)
11  SAVE /hiparnt/
12  common/ranseed/nseed
13  SAVE /ranseed/
14  common/histrng/nfp(300,15),pp(300,15),nft(300,15),pt(300,15)
15  SAVE /histrng/
16  IF(jp.EQ.0 .OR. jt.EQ.0) go to 25
17  DO 10 i=1,5
18  psc1(i)=pp(jp,i)
19  psc2(i)=pt(jt,i)
20 10 CONTINUE
21  CALL hijels(psc1,psc2)
22  dpp1=psc1(1)-pp(jp,1)
23  dpp2=psc1(2)-pp(jp,2)
24  dpt1=psc2(1)-pt(jt,1)
25  dpt2=psc2(2)-pt(jt,2)
26  pp(jp,6)=pp(jp,6)+dpp1/2.0
27  pp(jp,7)=pp(jp,7)+dpp2/2.0
28  pp(jp,8)=pp(jp,8)+dpp1/2.0
29  pp(jp,9)=pp(jp,9)+dpp2/2.0
30  pt(jt,6)=pt(jt,6)+dpt1/2.0
31  pt(jt,7)=pt(jt,7)+dpt2/2.0
32  pt(jt,8)=pt(jt,8)+dpt1/2.0
33  pt(jt,9)=pt(jt,9)+dpt2/2.0
34  DO 20 i=1,4
35  pp(jp,i)=psc1(i)
36  pt(jt,i)=psc2(i)
37 20 CONTINUE
38  nfp(jp,5)=max(1,nfp(jp,5))
39  nft(jt,5)=max(1,nft(jt,5))
40 C ********Perform elastic scattering between JP and JT
41  RETURN
42 C ********The following is for possible elastic cascade
43 c
44 25 IF(jp.EQ.0) go to 45
45  pabs=sqrt(pp(jp,1)**2+pp(jp,2)**2+pp(jp,3)**2)
46  bx=pp(jp,1)/pabs
47  by=pp(jp,2)/pabs
48  bz=pp(jp,3)/pabs
49  DO 40 i=1,ihnt2(1)
50  IF(i.EQ.jp) go to 40
51  dx=yp(1,i)-yp(1,jp)
52  dy=yp(2,i)-yp(2,jp)
53  dz=yp(3,i)-yp(3,jp)
54  dis=dx*bx+dy*by+dz*bz
55  IF(dis.LE.0) go to 40
56  bb=dx**2+dy**2+dz**2-dis**2
57  r2=bb*hipr1(40)/hipr1(31)/0.1
58 C ********mb=0.1*fm, YP is in fm,HIPR1(31) is in mb
59  gs=1.0-exp(-(hipr1(30)+hint1(11))/hipr1(31)/2.0
60  & *romg(r2))**2
61  gs0=1.0-exp(-(hipr1(30)+hint1(11))/hipr1(31)/2.0
62  & *romg(0.0))**2
63  IF(atl_ran(nseed).GT.gs/gs0) go to 40
64  DO 30 k=1,5
65  psc1(k)=pp(jp,k)
66  psc2(k)=pp(i,k)
67 30 CONTINUE
68  CALL hijels(psc1,psc2)
69  dpp1=psc1(1)-pp(jp,1)
70  dpp2=psc1(2)-pp(jp,2)
71  dpt1=psc2(1)-pp(i,1)
72  dpt2=psc2(2)-pp(i,2)
73  pp(jp,6)=pp(jp,6)+dpp1/2.0
74  pp(jp,7)=pp(jp,7)+dpp2/2.0
75  pp(jp,8)=pp(jp,8)+dpp1/2.0
76  pp(jp,9)=pp(jp,9)+dpp2/2.0
77  pp(i,6)=pp(i,6)+dpt1/2.0
78  pp(i,7)=pp(i,7)+dpt2/2.0
79  pp(i,8)=pp(i,8)+dpt1/2.0
80  pp(i,9)=pp(i,9)+dpt2/2.0
81  DO 35 k=1,5
82  pp(jp,k)=psc1(k)
83  pp(i,k)=psc2(k)
84 35 CONTINUE
85  nfp(i,5)=max(1,nfp(i,5))
86  go to 45
87 40 CONTINUE
88 45 IF(jt.EQ.0) go to 80
89 50 pabs=sqrt(pt(jt,1)**2+pt(jt,2)**2+pt(jt,3)**2)
90  bx=pt(jt,1)/pabs
91  by=pt(jt,2)/pabs
92  bz=pt(jt,3)/pabs
93  DO 70 i=1,ihnt2(3)
94  IF(i.EQ.jt) go to 70
95  dx=yt(1,i)-yt(1,jt)
96  dy=yt(2,i)-yt(2,jt)
97  dz=yt(3,i)-yt(3,jt)
98  dis=dx*bx+dy*by+dz*bz
99  IF(dis.LE.0) go to 70
100  bb=dx**2+dy**2+dz**2-dis**2
101  r2=bb*hipr1(40)/hipr1(31)/0.1
102 C ********mb=0.1*fm, YP is in fm,HIPR1(31) is in mb
103  gs=(1.0-exp(-(hipr1(30)+hint1(11))/hipr1(31)/2.0
104  & *romg(r2)))**2
105  gs0=(1.0-exp(-(hipr1(30)+hint1(11))/hipr1(31)/2.0
106  & *romg(0.0)))**2
107  IF(atl_ran(nseed).GT.gs/gs0) go to 70
108  DO 60 k=1,5
109  psc1(k)=pt(jt,k)
110  psc2(k)=pt(i,k)
111 60 CONTINUE
112  CALL hijels(psc1,psc2)
113  dpp1=psc1(1)-pt(jt,1)
114  dpp2=psc1(2)-pt(jt,2)
115  dpt1=psc2(1)-pt(i,1)
116  dpt2=psc2(2)-pt(i,2)
117  pt(jt,6)=pt(jt,6)+dpp1/2.0
118  pt(jt,7)=pt(jt,7)+dpp2/2.0
119  pt(jt,8)=pt(jt,8)+dpp1/2.0
120  pt(jt,9)=pt(jt,9)+dpp2/2.0
121  pt(i,6)=pt(i,6)+dpt1/2.0
122  pt(i,7)=pt(i,7)+dpt2/2.0
123  pt(i,8)=pt(i,8)+dpt1/2.0
124  pt(i,9)=pt(i,9)+dpt2/2.0
125  DO 65 k=1,5
126  pt(jt,k)=psc1(k)
127  pt(i,k)=psc2(k)
128 65 CONTINUE
129  nft(i,5)=max(1,nft(i,5))
130  go to 80
131 70 CONTINUE
132 80 RETURN
133  END