120a9922a93403acea584506db1d9fa653f7eafa
[u/mrichter/AliRoot.git] / HIJING / hijing1_36 / hijsft.F
1 * $Id$
2 C
3 C       
4 C*******************************************************************
5 C                                                                  *
6 C               Subroutine HIJSFT                                  *
7 C                                                                  *
8 C  Scatter two excited strings, JP from proj and JT from target    *
9 C*******************************************************************
10         SUBROUTINE HIJSFT(JP,JT,JOUT,IERROR)
11 #include "hijcrdn.inc"
12 #include "hiparnt.inc"
13 #include "hijdat.inc"
14 #include "hijjet1.inc"
15 #include "hijjet2.inc"
16 #include "histrng.inc"
17 #include "dpmcom1.inc"
18 #include "dpmcom2.inc"
19         SAVE
20 C*******************************************************************
21 C       JOUT-> the number
22 C       of hard scatterings preceding this soft collision. 
23 C       IHNT2(13)-> 1=
24 C       double diffrac 2=single diffrac, 3=non-single diffrac.
25 C*******************************************************************
26         IERROR=0
27         JJP=JP
28         JJT=JT
29         NDPM=0
30         IOPMAIN=0
31         IF(JP.GT.IHNT2(1) .OR. JT.GT.IHNT2(3)) RETURN
32
33         EPP=PP(JP,4)+PP(JP,3)
34         EPM=PP(JP,4)-PP(JP,3)
35         ETP=PT(JT,4)+PT(JT,3)
36         ETM=PT(JT,4)-PT(JT,3)
37
38         WP=EPP+ETP
39         WM=EPM+ETM
40         SW=WP*WM
41 C               ********total W+,W- and center-of-mass energy
42
43         IF(WP.LT.0.0 .OR. WM.LT.0.0) GO TO 1000
44
45         IF(JOUT.EQ.0) THEN
46                 IF(EPP.LT.0.0) GO TO 1000
47                 IF(EPM.LT.0.0) GO TO 1000
48                 IF(ETP.LT.0.0) GO TO 1000
49                 IF(ETM.LT.0.0) GO TO 1000    
50                 IF(EPP/(EPM+0.01).LE.ETP/(ETM+0.01)) RETURN
51         ENDIF
52 C               ********For strings which does not follow a jet-prod,
53 C                       scatter only if Ycm(JP)>Ycm(JT). When jets
54 C                       are produced just before this collision
55 C                       this requirement has already be enforced
56 C                       (see SUBROUTINE HIJHRD)
57         IHNT2(11)=JP
58         IHNT2(12)=JT
59 C
60 C
61 C
62         MISS=0
63         PKC1=0.0
64         PKC2=0.0
65         PKC11=0.0
66         PKC12=0.0
67         PKC21=0.0
68         PKC22=0.0
69         DPKC11=0.0
70         DPKC12=0.0
71         DPKC21=0.0
72         DPKC22=0.0
73         IF(NFP(JP,10).EQ.1.OR.NFT(JT,10).EQ.1) THEN
74            IF(NFP(JP,10).EQ.1) THEN
75               PHI1=ULANGL_HIJING(PP(JP,10),PP(JP,11))
76               PPJET=SQRT(PP(JP,10)**2+PP(JP,11)**2)
77               PKC1=PPJET
78               PKC11=PP(JP,10)
79               PKC12=PP(JP,11)
80            ENDIF
81            IF(NFT(JT,10).EQ.1) THEN
82               PHI2=ULANGL_HIJING(PT(JT,10),PT(JT,11))
83               PTJET=SQRT(PT(JT,10)**2+PT(JT,11)**2)
84               PKC2=PTJET
85               PKC21=PT(JT,10)
86               PKC22=PT(JT,11)
87            ENDIF
88            IF(IHPR2(4).GT.0.AND.IHNT2(1).GT.1.AND.IHNT2(3).GT.1) THEN
89               IF(NFP(JP,10).EQ.0) THEN
90                  PHI=-PHI2
91               ELSE IF(NFT(JT,10).EQ.0) THEN
92                  PHI=PHI1
93               ELSE
94                  PHI=(PHI1+PHI2-HIPR1(40))/2.0
95               ENDIF
96               BX=HINT1(19)*COS(HINT1(20))
97               BY=HINT1(19)*SIN(HINT1(20))
98               XP0=YP(1,JP)
99               YP0=YP(2,JP)
100               XT0=YT(1,JT)+BX
101               YT0=YT(2,JT)+BY
102               R1=MAX(1.2*IHNT2(1)**0.3333333,
103      &               SQRT(XP0**2+YP0**2))
104               R2=MAX(1.2*IHNT2(3)**0.3333333,
105      &               SQRT((XT0-BX)**2+(YT0-BY)**2))
106               IF(ABS(COS(PHI)).LT.1.0E-5) THEN
107                  DD1=R1
108                  DD2=R1
109 C PH             DD3=ABS(BY+SQRT(R2**2-(XP0-BX)**2)-YP0)
110 C PH             DD4=ABS(BY-SQRT(R2**2-(XP0-BX)**2)-YP0)
111                  DD3=ABS(BY+SQRT(R2**2-(XT0-BX)**2)-YP0)
112                  DD4=ABS(BY-SQRT(R2**2-(XT0-BX)**2)-YP0)
113                  GO TO 5
114               ENDIF
115               BB=2.0*SIN(PHI)*(COS(PHI)*YP0-SIN(PHI)*XP0)
116               CC=(YP0**2-R1**2)*COS(PHI)**2+XP0*SIN(PHI)*(
117      &                          XP0*SIN(PHI)-2.0*YP0*COS(PHI))
118               DD=BB**2-4.0*CC
119               IF(DD.LT.0.0) GO TO 10
120               XX1=(-BB+SQRT(DD))/2.0
121               XX2=(-BB-SQRT(DD))/2.0
122               DD1=ABS((XX1-XP0)/COS(PHI))
123               DD2=ABS((XX2-XP0)/COS(PHI))
124 C                       
125               BB=2.0*SIN(PHI)*(COS(PHI)*(YT0-BY)-SIN(PHI)*XT0)-2.0*BX
126               CC=(BX**2+(YT0-BY)**2-R2**2)*COS(PHI)**2+XT0*SIN(PHI)
127      &           *(XT0*SIN(PHI)-2.0*COS(PHI)*(YT0-BY))
128      &           -2.0*BX*SIN(PHI)*(COS(PHI)*(YT0-BY)-SIN(PHI)*XT0)
129               DD=BB**2-4.0*CC
130               IF(DD.LT.0.0) GO TO 10
131               XX1=(-BB+SQRT(DD))/2.0
132               XX2=(-BB-SQRT(DD))/2.0
133               DD3=ABS((XX1-XT0)/COS(PHI))
134               DD4=ABS((XX2-XT0)/COS(PHI))
135 C
136  5            DD1=MIN(DD1,DD3)
137               DD2=MIN(DD2,DD4)
138               IF(DD1.LT.HIPR1(13)) DD1=0.0
139               IF(DD2.LT.HIPR1(13)) DD2=0.0
140               IF(NFP(JP,10).EQ.1.AND.PPJET.GT.HIPR1(11)) THEN
141                  IF (IHPR2(50) .EQ. 1) THEN 
142                     DEDX0 = HIPR1(14)*LOG10(PPJET)/LOG10(5.)
143                  ELSE
144                     DEDX0 = HIPR1(14)
145                  ENDIF
146                  DP1=DD1*DEDX0/2.0
147                  DP1=MIN(DP1,PPJET-HIPR1(11))
148                  PKC1=PPJET-DP1
149                  DPX1=COS(PHI1)*DP1
150                  DPY1=SIN(PHI1)*DP1
151                  PKC11=PP(JP,10)-DPX1
152                  PKC12=PP(JP,11)-DPY1
153                  IF(DP1.GT.0.0) THEN
154                     CTHEP=PP(JP,12)/SQRT(PP(JP,12)**2+PPJET**2)
155                     DPZ1=DP1*CTHEP/SQRT(1.0-CTHEP**2)
156                     DPE1=SQRT(DPX1**2+DPY1**2+DPZ1**2)
157                     EPPPRM=PP(JP,4)+PP(JP,3)-DPE1-DPZ1
158                     EPMPRM=PP(JP,4)-PP(JP,3)-DPE1+DPZ1
159                     IF(EPPPRM.LE.0.0.OR.EPMPRM.LE.0.0) GO TO 15
160                     EPP=EPPPRM
161                     EPM=EPMPRM
162                     PP(JP,10)=PKC11
163                     PP(JP,11)=PKC12
164                     NPJ(JP)=NPJ(JP)+1
165                     KFPJ(JP,NPJ(JP))=21
166                     PJPX(JP,NPJ(JP))=DPX1
167                     PJPY(JP,NPJ(JP))=DPY1
168                     PJPZ(JP,NPJ(JP))=DPZ1
169                     PJPE(JP,NPJ(JP))=DPE1
170                     PJPM(JP,NPJ(JP))=0.0
171                     PP(JP,3)=PP(JP,3)-DPZ1
172                     PP(JP,4)=PP(JP,4)-DPE1
173                  ENDIF
174               ENDIF
175  15           IF(NFT(JT,10).EQ.1.AND.PTJET.GT.HIPR1(11)) THEN
176                  IF (IHPR2(50) .EQ. 1) THEN 
177                     DEDX0 = HIPR1(14)*LOG10(PTJET)/LOG10(5.)
178                  ELSE
179                     DEDX0 = HIPR1(14)
180                  ENDIF
181                  DP2=DD2*DEDX0/2.0
182                  DP2=MIN(DP2,PTJET-HIPR1(11))
183                  PKC2=PTJET-DP2
184                  DPX2=COS(PHI2)*DP2
185                  DPY2=SIN(PHI2)*DP2
186                  PKC21=PT(JT,10)-DPX2
187                  PKC22=PT(JT,11)-DPY2
188                  IF(DP2.GT.0.0) THEN
189                     CTHET=PT(JT,12)/SQRT(PT(JT,12)**2+PTJET**2)
190                     DPZ2=DP2*CTHET/SQRT(1.0-CTHET**2)
191                     DPE2=SQRT(DPX2**2+DPY2**2+DPZ2**2)
192                     ETPPRM=PT(JT,4)+PT(JT,3)-DPE2-DPZ2
193                     ETMPRM=PT(JT,4)-PT(JT,3)-DPE2+DPZ2
194                     IF(ETPPRM.LE.0.0.OR.ETMPRM.LE.0.0) GO TO 16
195                     ETP=ETPPRM
196                     ETM=ETMPRM
197                     PT(JT,10)=PKC21
198                     PT(JT,11)=PKC22
199                     NTJ(JT)=NTJ(JT)+1
200                     KFTJ(JT,NTJ(JT))=21
201                     PJTX(JT,NTJ(JT))=DPX2
202                     PJTY(JT,NTJ(JT))=DPY2
203                     PJTZ(JT,NTJ(JT))=DPZ2
204                     PJTE(JT,NTJ(JT))=DPE2
205                     PJTM(JT,NTJ(JT))=0.0
206                     PT(JT,3)=PT(JT,3)-DPZ2
207                     PT(JT,4)=PT(JT,4)-DPE2
208                  ENDIF
209               ENDIF
210  16           DPKC11=-(PP(JP,10)-PKC11)/2.0
211               DPKC12=-(PP(JP,11)-PKC12)/2.0
212               DPKC21=-(PT(JT,10)-PKC21)/2.0
213               DPKC22=-(PT(JT,11)-PKC22)/2.0
214               WP=EPP+ETP
215               WM=EPM+ETM
216               SW=WP*WM
217            ENDIF
218         ENDIF
219 C               ********If jet is quenched the pt from valence quark
220 C                       hard scattering has to reduced by d*kapa
221 C
222 C   
223 10      PTP02=PP(JP,1)**2+PP(JP,2)**2
224         PTT02=PT(JT,1)**2+PT(JT,2)**2
225 C       
226         AMQ=MAX(PP(JP,14)+PP(JP,15),PT(JT,14)+PT(JT,15))
227         AMX=HIPR1(1)+AMQ
228 C               ********consider mass cut-off for strings which
229 C                       must also include quark's mass
230         AMP0=AMX
231         DPM0=AMX
232         NFDP=0
233         IF(NFP(JP,5).LE.2.AND.NFP(JP,3).NE.0) THEN
234                 AMP0=ULMASS_HIJING(NFP(JP,3))
235                 NFDP=NFP(JP,3)+2*NFP(JP,3)/ABS(NFP(JP,3))
236                 DPM0=ULMASS_HIJING(NFDP)
237                 IF(DPM0.LE.0.0) THEN
238                         NFDP=NFDP-2*NFDP/ABS(NFDP)
239                         DPM0=ULMASS_HIJING(NFDP)
240                 ENDIF
241         ENDIF
242         AMT0=AMX
243         DTM0=AMX
244         NFDT=0
245         IF(NFT(JT,5).LE.2.AND.NFT(JT,3).NE.0) THEN
246                 AMT0=ULMASS_HIJING(NFT(JT,3))
247                 NFDT=NFT(JT,3)+2*NFT(JT,3)/ABS(NFT(JT,3))
248                 DTM0=ULMASS_HIJING(NFDT)
249                 IF(DTM0.LE.0.0) THEN
250                         NFDT=NFDT-2*NFDT/ABS(NFDT)
251                         DTM0=ULMASS_HIJING(NFDT)
252                 ENDIF
253         ENDIF
254 C       
255         AMPN=SQRT(AMP0**2+PTP02)
256         AMTN=SQRT(AMT0**2+PTT02)
257         SNN=(AMPN+AMTN)**2+0.001
258 C
259         IF(SW.LT.SNN+0.001) GO TO 4000
260 C               ********Scatter only if SW>SNN
261 C*****give some PT kick to the two exited strings******************
262 20      SWPTN=4.0*(MAX(AMP0,AMT0)**2+MAX(PTP02,PTT02))
263         SWPTD=4.0*(MAX(DPM0,DTM0)**2+MAX(PTP02,PTT02))
264         SWPTX=4.0*(AMX**2+MAX(PTP02,PTT02))
265         IF(SW.LE.SWPTN) THEN
266                 PKCMX=0.0
267         ELSE IF(SW.GT.SWPTN .AND. SW.LE.SWPTD
268      &          .AND.NPJ(JP).EQ.0.AND.NTJ(JT).EQ.0) THEN
269            PKCMX=SQRT(SW/4.0-MAX(AMP0,AMT0)**2)
270      &           -SQRT(MAX(PTP02,PTT02))
271         ELSE IF(SW.GT.SWPTD .AND. SW.LE.SWPTX
272      &          .AND.NPJ(JP).EQ.0.AND.NTJ(JT).EQ.0) THEN
273            PKCMX=SQRT(SW/4.0-MAX(DPM0,DTM0)**2)
274      &           -SQRT(MAX(PTP02,PTT02))
275         ELSE IF(SW.GT.SWPTX) THEN
276            PKCMX=SQRT(SW/4.0-AMX**2)-SQRT(MAX(PTP02,PTT02))
277         ENDIF
278 C               ********maximun PT kick
279 C*********************************************************
280 C
281         IF(NFP(JP,10).EQ.1.OR.NFT(JT,10).EQ.1) THEN
282                 IF(PKC1.GT.PKCMX) THEN
283                         PKC1=PKCMX
284                         PKC11=PKC1*COS(PHI1)
285                         PKC12=PKC1*SIN(PHI1)
286                         DPKC11=-(PP(JP,10)-PKC11)/2.0
287                         DPKC12=-(PP(JP,11)-PKC12)/2.0
288                 ENDIF
289                 IF(PKC2.GT.PKCMX) THEN
290                         PKC2=PKCMX
291                         PKC21=PKC2*COS(PHI2)
292                         PKC22=PKC2*SIN(PHI2)
293                         DPKC21=-(PT(JT,10)-PKC21)/2.0
294                         DPKC22=-(PT(JT,11)-PKC22)/2.0
295                 ENDIF
296                 DPKC1=DPKC11+DPKC21
297                 DPKC2=DPKC12+DPKC22
298                 NFP(JP,10)=-NFP(JP,10)
299                 NFT(JT,10)=-NFT(JT,10)
300                 GO TO 40
301         ENDIF
302 C               ********If the valence quarks had a hard-collision
303 C                       the pt kick is the pt from hard-collision.
304         I_SNG=0
305         IF(IHPR2(13).NE.0 .AND. RLU_HIJING(0).LE.HIDAT(4)) I_SNG=1
306         IF((NFP(JP,5).EQ.3 .OR.NFT(JT,5).EQ.3).OR.
307      &          (NPJ(JP).NE.0.OR.NFP(JP,10).NE.0).OR.
308      &          (NTJ(JT).NE.0.OR.NFT(JT,10).NE.0)) I_SNG=0
309 C
310 C               ********decite whether to have single-diffractive
311         IF(IHPR2(5).EQ.0) THEN
312                 PKC=HIPR1(2)*SQRT(-ALOG(1.0-RLU_HIJING(0)
313      &                  *(1.0-EXP(-PKCMX**2/HIPR1(2)**2))))
314                 GO TO 30
315         ENDIF
316         PKC=HIRND2(3,0.0,PKCMX**2)
317         PKC=SQRT(PKC)
318         IF(PKC.GT.HIPR1(20)) 
319      &     PKC=HIPR1(2)*SQRT(-ALOG(EXP(-HIPR1(20)**2/HIPR1(2)**2)
320      &         -RLU_HIJING(0)*(EXP(-HIPR1(20)**2/HIPR1(2)**2)-
321      &         EXP(-PKCMX**2/HIPR1(2)**2))))
322 C
323         IF(I_SNG.EQ.1) PKC=0.65*SQRT(
324      &          -ALOG(1.0-RLU_HIJING(0)*(1.0-EXP(-PKCMX**2/0.65**2))))
325 C                       ********select PT kick
326 30      PHI0=2.0*HIPR1(40)*RLU_HIJING(0)
327         PKC11=PKC*SIN(PHI0)
328         PKC12=PKC*COS(PHI0)
329         PKC21=-PKC11
330         PKC22=-PKC12
331         DPKC1=0.0
332         DPKC2=0.0
333 40      PP11=PP(JP,1)+PKC11-DPKC1
334         PP12=PP(JP,2)+PKC12-DPKC2
335         PT11=PT(JT,1)+PKC21-DPKC1
336         PT12=PT(JT,2)+PKC22-DPKC2
337         PTP2=PP11**2+PP12**2
338         PTT2=PT11**2+PT12**2
339 C
340         AMPN=SQRT(AMP0**2+PTP2)
341         AMTN=SQRT(AMT0**2+PTT2)
342         SNN=(AMPN+AMTN)**2+0.001
343 C***************************************
344         WP=EPP+ETP
345         WM=EPM+ETM
346         SW=WP*WM
347 C****************************************
348         IF(SW.LT.SNN) THEN
349            MISS=MISS+1
350            IF(MISS.LE.100) then
351               PKC=0.0
352               GO TO 30
353            ENDIF
354            IF(IHPR2(10).NE.0) 
355      &       WRITE(6,*) 'Error occured in Pt kick section of HIJSFT'
356            GO TO 4000
357         ENDIF
358 C******************************************************************
359         AMPD=SQRT(DPM0**2+PTP2)
360         AMTD=SQRT(DTM0**2+PTT2)
361
362         AMPX=SQRT(AMX**2+PTP2)
363         AMTX=SQRT(AMX**2+PTT2)
364
365         DPN=AMPN**2/SW
366         DTN=AMTN**2/SW
367         DPD=AMPD**2/SW
368         DTD=AMTD**2/SW
369         DPX=AMPX**2/SW
370         DTX=AMTX**2/SW
371 C
372         SPNTD=(AMPN+AMTD)**2
373         SPNTX=(AMPN+AMTX)**2
374 C                       ********CM energy if proj=N,targ=N*
375         SPDTN=(AMPD+AMTN)**2
376         SPXTN=(AMPX+AMTN)**2
377 C                       ********CM energy if proj=N*,targ=N
378         SPDTX=(AMPD+AMTX)**2
379         SPXTD=(AMPX+AMTD)**2
380         SDD=(AMPD+AMTD)**2
381         SXX=(AMPX+AMTX)**2
382
383 C
384 C       
385 C               ********CM energy if proj=delta, targ=delta
386 C****************There are many different cases**********
387 c       IF(IHPR2(15).EQ.1) GO TO 500
388 C
389 C               ********to have DPM type soft interactions
390 C
391  45     CONTINUE
392         IF(SW.GT.SXX+0.001) THEN
393            IF(I_SNG.EQ.0) THEN
394               D1=DPX
395               D2=DTX
396               NFP3=0
397               NFT3=0
398               GO TO 400
399            ELSE
400 c**** 5/30/1998 this is identical to the above statement. Added to
401 c**** avoid questional branching to block.
402               IF((NFP(JP,5).EQ.3 .AND.NFT(JT,5).EQ.3).OR.
403      &           (NPJ(JP).NE.0.OR.NFP(JP,10).NE.0).OR.
404      &           (NTJ(JT).NE.0.OR.NFT(JT,10).NE.0)) THEN
405                  D1=DPX
406                  D2=DTX
407                  NFP3=0
408                  NFT3=0
409                  GO TO 400
410               ENDIF
411 C               ********do not allow excited strings to have 
412 C                       single-diffr 
413               IF(RLU_HIJING(0).GT.0.5.OR.(NFT(JT,5).GT.2.OR.
414      &                NTJ(JT).NE.0.OR.NFT(JT,10).NE.0)) THEN
415                  D1=DPN
416                  D2=DTX
417                  NFP3=NFP(JP,3)
418                  NFT3=0
419                  GO TO 220
420               ELSE
421                  D1=DPX
422                  D2=DTN
423                  NFP3=0
424                  NFT3=NFT(JT,3)
425                  GO TO 240
426               ENDIF
427 C               ********have single diffractive collision
428            ENDIF
429         ELSE IF(SW.GT.MAX(SPDTX,SPXTD)+0.001 .AND.
430      &                          SW.LE.SXX+0.001) THEN
431            IF(((NPJ(JP).EQ.0.AND.NTJ(JT).EQ.0.AND.
432      &         RLU_HIJING(0).GT.0.5).OR.(NPJ(JP).EQ.0
433      &         .AND.NTJ(JT).NE.0)).AND.NFP(JP,5).LE.2) THEN
434               D1=DPD
435               D2=DTX
436               NFP3=NFDP
437               NFT3=0
438               GO TO 220
439            ELSE IF(NTJ(JT).EQ.0.AND.NFT(JT,5).LE.2) THEN
440               D1=DPX
441               D2=DTD
442               NFP3=0
443               NFT3=NFDT
444               GO TO 240
445            ENDIF
446            GO TO 4000
447         ELSE IF(SW.GT.MIN(SPDTX,SPXTD)+0.001.AND.
448      &                  SW.LE.MAX(SPDTX,SPXTD)+0.001) THEN
449            IF(SPDTX.LE.SPXTD.AND.NPJ(JP).EQ.0
450      &                       .AND.NFP(JP,5).LE.2) THEN
451               D1=DPD
452               D2=DTX
453               NFP3=NFDP
454               NFT3=0
455               GO TO 220
456            ELSE IF(SPDTX.GT.SPXTD.AND.NTJ(JT).EQ.0
457      &                       .AND.NFT(JT,5).LE.2) THEN
458               D1=DPX
459               D2=DTD
460               NFP3=0
461               NFT3=NFDT
462               GO TO 240
463            ENDIF
464 c*** 5/30/1998 added to avoid questional branching to another block
465 c*** this is identical to the statement following the next ELSE IF
466            IF(((NPJ(JP).EQ.0.AND.NTJ(JT).EQ.0
467      &       .AND.RLU_HIJING(0).GT.0.5).OR.(NPJ(JP).EQ.0
468      &        .AND.NTJ(JT).NE.0)).AND.NFP(JP,5).LE.2) THEN
469               D1=DPN
470               D2=DTX
471               NFP3=NFP(JP,3)
472               NFT3=0
473               GO TO 220
474            ELSE IF(NTJ(JT).EQ.0.AND.NFT(JT,5).LE.2) THEN
475               D1=DPX
476               D2=DTN
477               NFP3=0
478               NFT3=NFT(JT,3)
479               GO TO 240
480            ENDIF
481            GO TO 4000
482         ELSE IF(SW.GT.MAX(SPNTX,SPXTN)+0.001 .AND.
483      &                  SW.LE.MIN(SPDTX,SPXTD)+0.001) THEN
484            IF(((NPJ(JP).EQ.0.AND.NTJ(JT).EQ.0
485      &       .AND.RLU_HIJING(0).GT.0.5).OR.(NPJ(JP).EQ.0
486      &        .AND.NTJ(JT).NE.0)).AND.NFP(JP,5).LE.2) THEN
487               D1=DPN
488               D2=DTX
489               NFP3=NFP(JP,3)
490               NFT3=0
491               GO TO 220
492            ELSE IF(NTJ(JT).EQ.0.AND.NFT(JT,5).LE.2) THEN
493               D1=DPX
494               D2=DTN
495               NFP3=0
496               NFT3=NFT(JT,3)
497               GO TO 240
498            ENDIF
499            GO TO 4000
500         ELSE IF(SW.GT.MIN(SPNTX,SPXTN)+0.001 .AND.
501      &                  SW.LE.MAX(SPNTX,SPXTN)+0.001) THEN
502            IF(SPNTX.LE.SPXTN.AND.NPJ(JP).EQ.0
503      &                           .AND.NFP(JP,5).LE.2) THEN
504               D1=DPN
505               D2=DTX
506               NFP3=NFP(JP,3)
507               NFT3=0
508               GO TO 220
509            ELSEIF(SPNTX.GT.SPXTN.AND.NTJ(JT).EQ.0
510      &                           .AND.NFT(JT,5).LE.2) THEN
511               D1=DPX
512               D2=DTN
513               NFP3=0
514               NFT3=NFT(JT,3)
515               GO TO 240
516            ENDIF
517            GO TO 4000
518         ELSE IF(SW.LE.MIN(SPNTX,SPXTN)+0.001 .AND.
519      &                  (NPJ(JP).NE.0 .OR.NTJ(JT).NE.0)) THEN
520            GO TO 4000
521         ELSE IF(SW.LE.MIN(SPNTX,SPXTN)+0.001 .AND.
522      &          NFP(JP,5).GT.2.AND.NFT(JT,5).GT.2) THEN
523            GO TO 4000
524         ELSE IF(SW.GT.SDD+0.001.AND.SW.LE.
525      &                     MIN(SPNTX,SPXTN)+0.001) THEN
526            D1=DPD
527            D2=DTD
528            NFP3=NFDP
529            NFT3=NFDT
530            GO TO 100
531         ELSE IF(SW.GT.MAX(SPNTD,SPDTN)+0.001 
532      &                      .AND. SW.LE.SDD+0.001) THEN
533            IF(RLU_HIJING(0).GT.0.5) THEN
534               D1=DPD
535               D2=DTN
536               NFP3=NFDP
537               NFT3=NFT(JT,3)
538               GO TO 100
539            ELSE
540               D1=DPN
541               D2=DTD
542               NFP3=NFP(JP,3)
543               NFT3=NFDT
544               GO TO 100
545            ENDIF
546         ELSE IF(SW.GT.MIN(SPNTD,SPDTN)+0.001
547      &          .AND. SW.LE.MAX(SPNTD,SPDTN)+0.001) THEN
548            IF(SPNTD.GT.SPDTN) THEN
549               D1=DPD
550               D2=DTN
551               NFP3=NFDP
552               NFT3=NFT(JT,3)
553               GO TO 100
554            ELSE
555               D1=DPN
556               D2=DTD
557               NFP3=NFP(JP,3)
558               NFT3=NFDT
559               GO TO 100
560            ENDIF
561         ELSE IF(SW.LE.MIN(SPNTD,SPDTN)+0.001) THEN
562            D1=DPN
563            D2=DTN
564            NFP3=NFP(JP,3)
565            NFT3=NFT(JT,3)
566            GO TO 100
567         ENDIF
568         WRITE(6,*) ' Error in HIJSFT: There is no path to here'
569         RETURN
570 C
571 C***************  elastic scattering ***************
572 C       this is like elastic, both proj and targ mass
573 C       must be fixed
574 C***************************************************
575 100     NFP5=MAX(2,NFP(JP,5))
576         NFT5=MAX(2,NFT(JT,5))
577         BB1=1.0+D1-D2
578         BB2=1.0+D2-D1
579         IF(BB1**2.LT.4.0*D1 .OR. BB2**2.LT.4.0*D2) THEN
580                 MISS=MISS+1
581                 IF(MISS.GT.100.OR.PKC.EQ.0.0) GO TO 3000
582                 PKC=PKC*0.5
583                 GO TO 30
584         ENDIF
585         IF(RLU_HIJING(0).LT.0.5) THEN
586                 X1=(BB1-SQRT(BB1**2-4.0*D1))/2.0
587                 X2=(BB2-SQRT(BB2**2-4.0*D2))/2.0
588         ELSE
589                 X1=(BB1+SQRT(BB1**2-4.0*D1))/2.0
590                 X2=(BB2+SQRT(BB2**2-4.0*D2))/2.0
591         ENDIF
592         IHNT2(13)=2
593         GO TO 600
594 C
595 C********** Single diffractive ***********************
596 C either proj or targ's mass is fixed
597 C*****************************************************
598 220     NFP5=MAX(2,NFP(JP,5))
599         NFT5=3
600         IF(NFP3.EQ.0) NFP5=3
601         BB2=1.0+D2-D1
602         IF(BB2**2.LT.4.0*D2) THEN
603                 MISS=MISS+1
604                 IF(MISS.GT.100.OR.PKC.EQ.0.0) GO TO 3000
605                 PKC=PKC*0.5
606                 GO TO 30
607         ENDIF
608         XMIN=(BB2-SQRT(BB2**2-4.0*D2))/2.0
609         XMAX=(BB2+SQRT(BB2**2-4.0*D2))/2.0
610         MISS4=0
611 222     X2=HIRND2(6,XMIN,XMAX)
612         X1=D1/(1.0-X2)
613         IF(X2*(1.0-X1).LT.(D2+1.E-4/SW)) THEN
614                 MISS4=MISS4+1
615                 IF(MISS4.LE.1000) GO TO 222
616                 GO TO 5000
617         ENDIF
618         IHNT2(13)=2
619         GO TO 600
620 C                       ********Fix proj mass*********
621 240     NFP5=3
622         NFT5=MAX(2,NFT(JT,5))
623         IF(NFT3.EQ.0) NFT5=3
624         BB1=1.0+D1-D2
625         IF(BB1**2.LT.4.0*D1) THEN
626                 MISS=MISS+1
627                 IF(MISS.GT.100.OR.PKC.EQ.0.0) GO TO 3000
628                 PKC=PKC*0.5
629                 GO TO 30
630         ENDIF
631         XMIN=(BB1-SQRT(BB1**2-4.0*D1))/2.0
632         XMAX=(BB1+SQRT(BB1**2-4.0*D1))/2.0
633         MISS4=0
634 242     X1=HIRND2(6,XMIN,XMAX)
635         X2=D2/(1.0-X1)
636         IF(X1*(1.0-X2).LT.(D1+1.E-4/SW)) THEN
637                 MISS4=MISS4+1
638                 IF(MISS4.LE.1000) GO TO 242
639                 GO TO 5000
640         ENDIF
641         IHNT2(13)=2
642         GO TO 600
643 C                       ********Fix targ mass*********
644 C
645 C*************non-single diffractive**********************
646 C       both proj and targ may not be fixed in mass 
647 C*********************************************************
648 C
649 400     NFP5=3
650         NFT5=3
651         BB1=1.0+D1-D2
652         BB2=1.0+D2-D1
653         IF(BB1**2.LT.4.0*D1 .OR. BB2**2.LT.4.0*D2) THEN
654                 MISS=MISS+1
655                 IF(MISS.GT.100.OR.PKC.EQ.0.0) GO TO 3000
656                 PKC=PKC*0.5
657                 GO TO 30
658         ENDIF
659         XMIN1=(BB1-SQRT(BB1**2-4.0*D1))/2.0
660         XMAX1=(BB1+SQRT(BB1**2-4.0*D1))/2.0
661         XMIN2=(BB2-SQRT(BB2**2-4.0*D2))/2.0
662         XMAX2=(BB2+SQRT(BB2**2-4.0*D2))/2.0
663         MISS4=0 
664 410     X1=HIRND2(4,XMIN1,XMAX1)
665         X2=HIRND2(4,XMIN2,XMAX2)
666         IF(NFP(JP,5).EQ.3.OR.NFT(JT,5).EQ.3) THEN
667                 X1=HIRND2(6,XMIN1,XMAX1)
668                 X2=HIRND2(6,XMIN2,XMAX2)
669         ENDIF
670 C                       ********
671         IF(ABS(NFP(JP,1)*NFP(JP,2)).GT.1000000.OR.
672      &                  ABS(NFP(JP,1)*NFP(JP,2)).LT.100) THEN
673                 X1=HIRND2(5,XMIN1,XMAX1)
674         ENDIF
675         IF(ABS(NFT(JT,1)*NFT(JT,2)).GT.1000000.OR.
676      &                  ABS(NFT(JT,1)*NFT(JT,2)).LT.100) THEN
677                 X2=HIRND2(5,XMIN2,XMAX2)
678         ENDIF
679 c       IF(IOPMAIN.EQ.3) X1=HIRND2(6,XMIN1,XMAX1)
680 c       IF(IOPMAIN.EQ.2) X2=HIRND2(6,XMIN2,XMAX2) 
681 C       ********For q-qbar or (qq)-(qq)bar system use symetric
682 C               distribution, for q-(qq) or qbar-(qq)bar use
683 C               unsymetrical distribution
684 C
685         IF(ABS(NFP(JP,1)*NFP(JP,2)).GT.1000000) X1=1.0-X1
686         XXP=X1*(1.0-X2)
687         XXT=X2*(1.0-X1)
688         IF(XXP.LT.(D1+1.E-4/SW) .OR. XXT.LT.(D2+1.E-4/SW)) THEN
689                 MISS4=MISS4+1
690                 IF(MISS4.LE.1000) GO TO 410
691                 GO TO 5000
692         ENDIF
693         IHNT2(13)=3
694 C***************************************************
695 C***************************************************
696 600     CONTINUE
697         IF(X1*(1.0-X2).LT.(AMPN**2-1.E-4)/SW.OR.
698      &                  X2*(1.0-X1).LT.(AMTN**2-1.E-4)/SW) THEN
699                 MISS=MISS+1
700                 IF(MISS.GT.100.OR.PKC.EQ.0.0) GO TO 2000
701                 PKC=0.0
702                 GO TO 30
703         ENDIF
704 C
705         EPP=(1.0-X2)*WP
706         EPM=X1*WM
707         ETP=X2*WP
708         ETM=(1.0-X1)*WM
709         PP(JP,3)=(EPP-EPM)/2.0
710         PP(JP,4)=(EPP+EPM)/2.0
711         IF(EPP*EPM-PTP2.LT.0.0) GO TO 6000
712         PP(JP,5)=SQRT(EPP*EPM-PTP2)
713         NFP(JP,3)=NFP3
714         NFP(JP,5)=NFP5
715
716         PT(JT,3)=(ETP-ETM)/2.0
717         PT(JT,4)=(ETP+ETM)/2.0
718         IF(ETP*ETM-PTT2.LT.0.0) GO TO 6000
719         PT(JT,5)=SQRT(ETP*ETM-PTT2)
720         NFT(JT,3)=NFT3
721         NFT(JT,5)=NFT5
722 C*****recoil PT from hard-inter is shared by two end-partons 
723 C       so that pt=p1+p2
724         PP(JP,1)=PP11-PKC11
725         PP(JP,2)=PP12-PKC12
726
727         KICKDIP=1
728         KICKDIT=1
729         IF(ABS(NFP(JP,1)*NFP(JP,2)).GT.1000000.OR.
730      &                  ABS(NFP(JP,1)*NFP(JP,2)).LT.100) THEN
731                 KICKDIP=0
732         ENDIF
733         IF(ABS(NFT(JT,1)*NFT(JT,2)).GT.1000000.OR.
734      &                  ABS(NFT(JT,1)*NFT(JT,2)).LT.100) THEN
735                 KICKDIT=0
736         ENDIF
737         IF((KICKDIP.EQ.0.AND.RLU_HIJING(0).LT.0.5)
738      &     .OR.(KICKDIP.NE.0.AND.RLU_HIJING(0)
739      &     .LT.0.5/(1.0+(PKC11**2+PKC12**2)/HIPR1(22)**2))) THEN
740            PP(JP,6)=(PP(JP,1)-PP(JP,6)-PP(JP,8)-DPKC1)/2.0+PP(JP,6)
741            PP(JP,7)=(PP(JP,2)-PP(JP,7)-PP(JP,9)-DPKC2)/2.0+PP(JP,7)
742            PP(JP,8)=(PP(JP,1)-PP(JP,6)-PP(JP,8)-DPKC1)/2.0
743      &              +PP(JP,8)+PKC11
744            PP(JP,9)=(PP(JP,2)-PP(JP,7)-PP(JP,9)-DPKC2)/2.0
745      &              +PP(JP,9)+PKC12
746         ELSE
747            PP(JP,8)=(PP(JP,1)-PP(JP,6)-PP(JP,8)-DPKC1)/2.0+PP(JP,8)
748            PP(JP,9)=(PP(JP,2)-PP(JP,7)-PP(JP,9)-DPKC2)/2.0+PP(JP,9)
749            PP(JP,6)=(PP(JP,1)-PP(JP,6)-PP(JP,8)-DPKC1)/2.0
750      &              +PP(JP,6)+PKC11
751            PP(JP,7)=(PP(JP,2)-PP(JP,7)-PP(JP,9)-DPKC2)/2.0
752      &              +PP(JP,7)+PKC12
753         ENDIF
754         PP(JP,1)=PP(JP,6)+PP(JP,8)
755         PP(JP,2)=PP(JP,7)+PP(JP,9)
756 C                               ********pt kick for proj
757         PT(JT,1)=PT11-PKC21
758         PT(JT,2)=PT12-PKC22
759         IF((KICKDIT.EQ.0.AND.RLU_HIJING(0).LT.0.5)
760      &     .OR.(KICKDIT.NE.0.AND.RLU_HIJING(0)
761      &     .LT.0.5/(1.0+(PKC21**2+PKC22**2)/HIPR1(22)**2))) THEN
762            PT(JT,6)=(PT(JT,1)-PT(JT,6)-PT(JT,8)-DPKC1)/2.0+PT(JT,6)
763            PT(JT,7)=(PT(JT,2)-PT(JT,7)-PT(JT,9)-DPKC2)/2.0+PT(JT,7)
764            PT(JT,8)=(PT(JT,1)-PT(JT,6)-PT(JT,8)-DPKC1)/2.0
765      &              +PT(JT,8)+PKC21
766            PT(JT,9)=(PT(JT,2)-PT(JT,7)-PT(JT,9)-DPKC2)/2.0
767      &              +PT(JT,9)+PKC22
768         ELSE
769            PT(JT,8)=(PT(JT,1)-PT(JT,6)-PT(JT,8)-DPKC1)/2.0+PT(JT,8)
770            PT(JT,9)=(PT(JT,2)-PT(JT,7)-PT(JT,9)-DPKC2)/2.0+PT(JT,9)
771            PT(JT,6)=(PT(JT,1)-PT(JT,6)-PT(JT,8)-DPKC1)/2.0
772      &              +PT(JT,6)+PKC21
773            PT(JT,7)=(PT(JT,2)-PT(JT,7)-PT(JT,9)-DPKC2)/2.0
774      &              +PT(JT,7)+PKC22
775         ENDIF
776         PT(JT,1)=PT(JT,6)+PT(JT,8)
777         PT(JT,2)=PT(JT,7)+PT(JT,9)
778 C                       ********pt kick for targ
779
780         IF(NPJ(JP).NE.0) NFP(JP,5)=3
781         IF(NTJ(JT).NE.0) NFT(JT,5)=3
782 C                       ********jets must be connected to string
783         IF(EPP/(EPM+0.0001).LT.ETP/(ETM+0.0001).AND.
784      &                  ABS(NFP(JP,1)*NFP(JP,2)).LT.1000000)THEN
785                 DO 620 JSB=1,15
786                 PSB=PP(JP,JSB)
787                 PP(JP,JSB)=PT(JT,JSB)
788                 PT(JT,JSB)=PSB
789                 NSB=NFP(JP,JSB)
790                 NFP(JP,JSB)=NFT(JT,JSB)
791                 NFT(JT,JSB)=NSB
792 620             CONTINUE
793 C               ********when Ycm(JP)<Ycm(JT) after the collision
794 C                       exchange the positions of the two   
795         ENDIF
796 C
797         RETURN
798 C**************************************************
799 C**************************************************
800 1000    IERROR=1
801         IF(IHPR2(10).EQ.0) RETURN
802         WRITE(6,*) '    Fatal HIJSFT start error,abandon this event'
803         WRITE(6,*) '    PROJ E+,E-,W+',EPP,EPM,WP
804         WRITE(6,*) '    TARG E+,E-,W-',ETP,ETM,WM
805         WRITE(6,*) '    W+*W-, (APN+ATN)^2',SW,SNN
806         RETURN
807 2000    IERROR=0
808         IF(IHPR2(10).EQ.0) RETURN
809         WRITE(6,*) '    (2)energy partition fail,'
810         WRITE(6,*) '    HIJSFT not performed, but continue'
811         WRITE(6,*) '    MP1,MPN',X1*(1.0-X2)*SW,AMPN**2
812         WRITE(6,*) '    MT2,MTN',X2*(1.0-X1)*SW,AMTN**2
813         RETURN
814 3000    IERROR=0
815         IF(IHPR2(10).EQ.0) RETURN
816         WRITE(6,*) '    (3)something is wrong with the pt kick, '
817         WRITE(6,*) '    HIJSFT not performed, but continue'
818         WRITE(6,*) '    D1=',D1,' D2=',D2,' SW=',SW
819         WRITE(6,*) '    HISTORY NFP5=',NFP(JP,5),' NFT5=',NFT(JT,5)
820         WRITE(6,*) '    THIS COLLISON NFP5=',NFP5, ' NFT5=',NFT5
821         WRITE(6,*) '    # OF JET IN PROJ',NPJ(JP),' IN TARG',NTJ(JT)
822         RETURN
823 4000    IERROR=0
824         IF(IHPR2(10).EQ.0) RETURN
825         WRITE(6,*) '    (4)unable to choose process, but not harmful'
826         WRITE(6,*) '    HIJSFT not performed, but continue'
827         WRITE(6,*) '    PTP=',SQRT(PTP2),' PTT=',SQRT(PTT2),' SW=',SW
828         WRITE(6,*) '    AMCUT=',AMX,' JP=',JP,' JT=',JT
829         WRITE(6,*) '    HISTORY NFP5=',NFP(JP,5),' NFT5=',NFT(JT,5)
830         RETURN
831 5000    IERROR=0
832         IF(IHPR2(10).EQ.0) RETURN
833         WRITE(6,*) '    energy partition failed(5),for limited try'
834         WRITE(6,*) '    HIJSFT not performed, but continue'
835         WRITE(6,*) '    NFP5=',NFP5,' NFT5=',NFT5
836         WRITE(6,*) '    D1',D1,' X1(1-X2)',X1*(1.0-X2)
837         WRITE(6,*) '    D2',D2,' X2(1-X1)',X2*(1.0-X1)
838         RETURN
839 6000    PKC=0.0
840         MISS=MISS+1
841         IF(MISS.LT.100) GO TO 30
842         IERROR=1
843         IF(IHPR2(10).EQ.0) RETURN
844         WRITE(6,*) ' ERROR OCCURED, HIJSFT NOT PERFORMED'
845         WRITE(6,*) ' Abort this event'
846         WRITE(6,*) 'MTP,PTP2',EPP*EPM,PTP2,'  MTT,PTT2',ETP*ETM,PTT2 
847         RETURN
848         END