[u/mrichter/AliRoot.git] / HIJING / hijing1_36 / hijcsc.F

* $Id$
C
C*******************************************************************
C	This subroutine performs elastic scatterings and possible 
C	elastic cascading within their own nuclei
c*******************************************************************
	SUBROUTINE HIJCSC(JP,JT)
	DIMENSION PSC1(5),PSC2(5)
#define BLANKET_SAVE
#include "hijcrdn.inc"
#include "hiparnt.inc"
#include "histrng.inc"
	SAVE
	IF(JP.EQ.0 .OR. JT.EQ.0) GO TO 25
	DO 10 I=1,5
	PSC1(I)=PP(JP,I)
	PSC2(I)=PT(JT,I)
10	CONTINUE
	CALL HIJELS(PSC1,PSC2)
	DPP1=PSC1(1)-PP(JP,1)
	DPP2=PSC1(2)-PP(JP,2)
	DPT1=PSC2(1)-PT(JT,1)
	DPT2=PSC2(2)-PT(JT,2)
	PP(JP,6)=PP(JP,6)+DPP1/2.0
	PP(JP,7)=PP(JP,7)+DPP2/2.0
	PP(JP,8)=PP(JP,8)+DPP1/2.0
	PP(JP,9)=PP(JP,9)+DPP2/2.0
	PT(JT,6)=PT(JT,6)+DPT1/2.0
	PT(JT,7)=PT(JT,7)+DPT2/2.0
	PT(JT,8)=PT(JT,8)+DPT1/2.0
	PT(JT,9)=PT(JT,9)+DPT2/2.0
	DO 20 I=1,4
	PP(JP,I)=PSC1(I)
	PT(JT,I)=PSC2(I)
20	CONTINUE
	NFP(JP,5)=MAX(1,NFP(JP,5))
	NFT(JT,5)=MAX(1,NFT(JT,5))
C		********Perform elastic scattering between JP and JT
	RETURN
C		********The following is for possible elastic cascade
c
25	IF(JP.EQ.0) GO TO 45
	PABS=SQRT(PP(JP,1)**2+PP(JP,2)**2+PP(JP,3)**2)
	BX=PP(JP,1)/PABS
	BY=PP(JP,2)/PABS
	BZ=PP(JP,3)/PABS
	DO 40 I=1,IHNT2(1)
		IF(I.EQ.JP) GO TO 40
		DX=YP(1,I)-YP(1,JP)
		DY=YP(2,I)-YP(2,JP)
		DZ=YP(3,I)-YP(3,JP)
		DIS=DX*BX+DY*BY+DZ*BZ
		IF(DIS.LE.0) GO TO 40
		BB=DX**2+DY**2+DZ**2-DIS**2
		R2=BB*HIPR1(40)/HIPR1(31)/0.1
C		********mb=0.1*fm, YP is in fm,HIPR1(31) is in mb
		GS=1.0-EXP(-(HIPR1(30)+HINT1(11))/HIPR1(31)/2.0
     &			*ROMG(R2))**2
		GS0=1.0-EXP(-(HIPR1(30)+HINT1(11))/HIPR1(31)/2.0
     &			*ROMG(0.0))**2
		IF(RLU_HIJING(0).GT.GS/GS0) GO TO 40
		DO 30 K=1,5
			PSC1(K)=PP(JP,K)
			PSC2(K)=PP(I,K)
30		CONTINUE
		CALL HIJELS(PSC1,PSC2)
		DPP1=PSC1(1)-PP(JP,1)
		DPP2=PSC1(2)-PP(JP,2)
		DPT1=PSC2(1)-PP(I,1)
		DPT2=PSC2(2)-PP(I,2)
		PP(JP,6)=PP(JP,6)+DPP1/2.0
		PP(JP,7)=PP(JP,7)+DPP2/2.0
		PP(JP,8)=PP(JP,8)+DPP1/2.0
		PP(JP,9)=PP(JP,9)+DPP2/2.0
		PP(I,6)=PP(I,6)+DPT1/2.0
		PP(I,7)=PP(I,7)+DPT2/2.0
		PP(I,8)=PP(I,8)+DPT1/2.0
		PP(I,9)=PP(I,9)+DPT2/2.0
		DO 35 K=1,5
			PP(JP,K)=PSC1(K)
			PP(I,K)=PSC2(K)
35		CONTINUE
		NFP(I,5)=MAX(1,NFP(I,5))
		GO TO 45
40	CONTINUE
45	IF(JT.EQ.0) GO TO 80
50	PABS=SQRT(PT(JT,1)**2+PT(JT,2)**2+PT(JT,3)**2)
	BX=PT(JT,1)/PABS
	BY=PT(JT,2)/PABS
	BZ=PT(JT,3)/PABS
	DO 70 I=1,IHNT2(3)
		IF(I.EQ.JT) GO TO 70
		DX=YT(1,I)-YT(1,JT)
		DY=YT(2,I)-YT(2,JT)
		DZ=YT(3,I)-YT(3,JT)
		DIS=DX*BX+DY*BY+DZ*BZ
		IF(DIS.LE.0) GO TO 70
		BB=DX**2+DY**2+DZ**2-DIS**2
		R2=BB*HIPR1(40)/HIPR1(31)/0.1
C		********mb=0.1*fm, YP is in fm,HIPR1(31) is in mb
		GS=(1.0-EXP(-(HIPR1(30)+HINT1(11))/HIPR1(31)/2.0
     &			*ROMG(R2)))**2
		GS0=(1.0-EXP(-(HIPR1(30)+HINT1(11))/HIPR1(31)/2.0
     &			*ROMG(0.0)))**2
		IF(RLU_HIJING(0).GT.GS/GS0) GO TO 70
		DO 60 K=1,5
			PSC1(K)=PT(JT,K)
			PSC2(K)=PT(I,K)
60		CONTINUE
		CALL HIJELS(PSC1,PSC2)
		DPP1=PSC1(1)-PT(JT,1)
		DPP2=PSC1(2)-PT(JT,2)
		DPT1=PSC2(1)-PT(I,1)
		DPT2=PSC2(2)-PT(I,2)
		PT(JT,6)=PT(JT,6)+DPP1/2.0
		PT(JT,7)=PT(JT,7)+DPP2/2.0
		PT(JT,8)=PT(JT,8)+DPP1/2.0
		PT(JT,9)=PT(JT,9)+DPP2/2.0
		PT(I,6)=PT(I,6)+DPT1/2.0
		PT(I,7)=PT(I,7)+DPT2/2.0
		PT(I,8)=PT(I,8)+DPT1/2.0
		PT(I,9)=PT(I,9)+DPT2/2.0
		DO 65 K=1,5
			PT(JT,K)=PSC1(K)
			PT(I,K)=PSC2(K)
65		CONTINUE
		NFT(I,5)=MAX(1,NFT(I,5))
		GO TO 80
70	CONTINUE
80	RETURN
	END
Commit	Line	Data
e74335a4	1	* $Id$
	2	C
	3	C*******************************************************************
	4	C This subroutine performs elastic scatterings and possible
	5	C elastic cascading within their own nuclei
	6	c*******************************************************************
	7	SUBROUTINE HIJCSC(JP,JT)
	8	DIMENSION PSC1(5),PSC2(5)
bc676b8e	9	#define BLANKET_SAVE
e74335a4	10	#include "hijcrdn.inc"
	11	#include "hiparnt.inc"
	12	#include "histrng.inc"
	13	SAVE
	14	IF(JP.EQ.0 .OR. JT.EQ.0) GO TO 25
	15	DO 10 I=1,5
	16	PSC1(I)=PP(JP,I)
	17	PSC2(I)=PT(JT,I)
	18	10 CONTINUE
	19	CALL HIJELS(PSC1,PSC2)
	20	DPP1=PSC1(1)-PP(JP,1)
	21	DPP2=PSC1(2)-PP(JP,2)
	22	DPT1=PSC2(1)-PT(JT,1)
	23	DPT2=PSC2(2)-PT(JT,2)
	24	PP(JP,6)=PP(JP,6)+DPP1/2.0
	25	PP(JP,7)=PP(JP,7)+DPP2/2.0
	26	PP(JP,8)=PP(JP,8)+DPP1/2.0
	27	PP(JP,9)=PP(JP,9)+DPP2/2.0
	28	PT(JT,6)=PT(JT,6)+DPT1/2.0
	29	PT(JT,7)=PT(JT,7)+DPT2/2.0
	30	PT(JT,8)=PT(JT,8)+DPT1/2.0
	31	PT(JT,9)=PT(JT,9)+DPT2/2.0
	32	DO 20 I=1,4
	33	PP(JP,I)=PSC1(I)
	34	PT(JT,I)=PSC2(I)
	35	20 CONTINUE
	36	NFP(JP,5)=MAX(1,NFP(JP,5))
	37	NFT(JT,5)=MAX(1,NFT(JT,5))
	38	C ********Perform elastic scattering between JP and JT
	39	RETURN
	40	C ********The following is for possible elastic cascade
	41	c
	42	25 IF(JP.EQ.0) GO TO 45
	43	PABS=SQRT(PP(JP,1)2+PP(JP,2)2+PP(JP,3)**2)
	44	BX=PP(JP,1)/PABS
	45	BY=PP(JP,2)/PABS
	46	BZ=PP(JP,3)/PABS
	47	DO 40 I=1,IHNT2(1)
	48	IF(I.EQ.JP) GO TO 40
	49	DX=YP(1,I)-YP(1,JP)
	50	DY=YP(2,I)-YP(2,JP)
	51	DZ=YP(3,I)-YP(3,JP)
	52	DIS=DXBX+DYBY+DZ*BZ
	53	IF(DIS.LE.0) GO TO 40
	54	BB=DX2+DY2+DZ2-DIS2
	55	R2=BB*HIPR1(40)/HIPR1(31)/0.1
	56	C *******mb=0.1fm, YP is in fm,HIPR1(31) is in mb
	57	GS=1.0-EXP(-(HIPR1(30)+HINT1(11))/HIPR1(31)/2.0
	58	& ROMG(R2))*2
	59	GS0=1.0-EXP(-(HIPR1(30)+HINT1(11))/HIPR1(31)/2.0
	60	& ROMG(0.0))*2
	61	IF(RLU_HIJING(0).GT.GS/GS0) GO TO 40
	62	DO 30 K=1,5
	63	PSC1(K)=PP(JP,K)
	64	PSC2(K)=PP(I,K)
	65	30 CONTINUE
	66	CALL HIJELS(PSC1,PSC2)
	67	DPP1=PSC1(1)-PP(JP,1)
	68	DPP2=PSC1(2)-PP(JP,2)
	69	DPT1=PSC2(1)-PP(I,1)
	70	DPT2=PSC2(2)-PP(I,2)
	71	PP(JP,6)=PP(JP,6)+DPP1/2.0
	72	PP(JP,7)=PP(JP,7)+DPP2/2.0
	73	PP(JP,8)=PP(JP,8)+DPP1/2.0
74	PP(JP,9)=PP(JP,9)+DPP2/2.0
75	PP(I,6)=PP(I,6)+DPT1/2.0
76	PP(I,7)=PP(I,7)+DPT2/2.0
77	PP(I,8)=PP(I,8)+DPT1/2.0
78	PP(I,9)=PP(I,9)+DPT2/2.0
79	DO 35 K=1,5
80	PP(JP,K)=PSC1(K)
81	PP(I,K)=PSC2(K)
82	35 CONTINUE
83	NFP(I,5)=MAX(1,NFP(I,5))
84	GO TO 45
85	40 CONTINUE
86	45 IF(JT.EQ.0) GO TO 80
87	50 PABS=SQRT(PT(JT,1)2+PT(JT,2)2+PT(JT,3)**2)
88	BX=PT(JT,1)/PABS
89	BY=PT(JT,2)/PABS
90	BZ=PT(JT,3)/PABS
91	DO 70 I=1,IHNT2(3)
92	IF(I.EQ.JT) GO TO 70
93	DX=YT(1,I)-YT(1,JT)
94	DY=YT(2,I)-YT(2,JT)
95	DZ=YT(3,I)-YT(3,JT)
96	DIS=DXBX+DYBY+DZ*BZ
97	IF(DIS.LE.0) GO TO 70
98	BB=DX2+DY2+DZ2-DIS2
99	R2=BB*HIPR1(40)/HIPR1(31)/0.1
100	C *******mb=0.1fm, YP is in fm,HIPR1(31) is in mb
101	GS=(1.0-EXP(-(HIPR1(30)+HINT1(11))/HIPR1(31)/2.0
102	& ROMG(R2)))*2
103	GS0=(1.0-EXP(-(HIPR1(30)+HINT1(11))/HIPR1(31)/2.0
104	& ROMG(0.0)))*2
105	IF(RLU_HIJING(0).GT.GS/GS0) GO TO 70
106	DO 60 K=1,5
107	PSC1(K)=PT(JT,K)
108	PSC2(K)=PT(I,K)
109	60 CONTINUE
110	CALL HIJELS(PSC1,PSC2)
111	DPP1=PSC1(1)-PT(JT,1)
112	DPP2=PSC1(2)-PT(JT,2)
113	DPT1=PSC2(1)-PT(I,1)
114	DPT2=PSC2(2)-PT(I,2)
115	PT(JT,6)=PT(JT,6)+DPP1/2.0
116	PT(JT,7)=PT(JT,7)+DPP2/2.0
117	PT(JT,8)=PT(JT,8)+DPP1/2.0
118	PT(JT,9)=PT(JT,9)+DPP2/2.0
119	PT(I,6)=PT(I,6)+DPT1/2.0
120	PT(I,7)=PT(I,7)+DPT2/2.0
121	PT(I,8)=PT(I,8)+DPT1/2.0
122	PT(I,9)=PT(I,9)+DPT2/2.0
123	DO 65 K=1,5
124	PT(JT,K)=PSC1(K)
125	PT(I,K)=PSC2(K)
126	65 CONTINUE
127	NFT(I,5)=MAX(1,NFT(I,5))
128	GO TO 80
129	70 CONTINUE
130	80 RETURN
131	END