Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
hijsft.f
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file hijsft.f
1 C
2 C
3 C*******************************************************************
4 C *
5 C Subroutine HIJSFT *
6 C *
7 C Scatter two excited strings, JP from proj and JT from target *
8 C*******************************************************************
9  SUBROUTINE hijsft(JP,JT,JOUT,IERROR)
10  common/hijcrdn/yp(3,300),yt(3,300)
11  SAVE /hijcrdn/
12  common/hiparnt/hipr1(100),ihpr2(50),hint1(100),ihnt2(50)
13  SAVE /hiparnt/
14  common/hijdat/hidat0(10,10),hidat(10)
15  SAVE /hijdat/
16  common/ranseed/nseed
17  SAVE /ranseed/
18  common/hijjet1/npj(300),kfpj(300,500),pjpx(300,500),
19  & pjpy(300,500),pjpz(300,500),pjpe(300,500),
20  & pjpm(300,500),ntj(300),kftj(300,500),
21  & pjtx(300,500),pjty(300,500),pjtz(300,500),
22  & pjte(300,500),pjtm(300,500)
23  SAVE /hijjet1/
24  common/hijjet2/nsg,njsg(900),iasg(900,3),k1sg(900,100),
25  & k2sg(900,100),pxsg(900,100),pysg(900,100),
26  & pzsg(900,100),pesg(900,100),pmsg(900,100)
27  SAVE /hijjet2/
28  common/histrng/nfp(300,15),pp(300,15),nft(300,15),pt(300,15)
29  SAVE /histrng/
30  common/dpmcom1/jjp,jjt,amp,amt,apx0,atx0,ampn,amtn,amp0,amt0,
31  & nfdp,nfdt,wp,wm,sw,xremp,xremt,dpkc1,dpkc2,pp11,pp12,
32  & pt11,pt12,ptp2,ptt2
33  SAVE /dpmcom1/
34  common/dpmcom2/ndpm,kdpm(20,2),pdpm1(20,5),pdpm2(20,5)
35  SAVE /dpmcom2/
36 C*******************************************************************
37 C JOUT-> the number
38 C of hard scatterings preceding this soft collision.
39 C IHNT2(13)-> 1=
40 C double diffrac 2=single diffrac, 3=non-single diffrac.
41 C*******************************************************************
42  ierror=0
43  jjp=jp
44  jjt=jt
45  ndpm=0
46  iopmain=0
47  IF(jp.GT.ihnt2(1) .OR. jt.GT.ihnt2(3)) RETURN
48 
49  epp=pp(jp,4)+pp(jp,3)
50  epm=pp(jp,4)-pp(jp,3)
51  etp=pt(jt,4)+pt(jt,3)
52  etm=pt(jt,4)-pt(jt,3)
53 
54  wp=epp+etp
55  wm=epm+etm
56  sw=wp*wm
57 C ********total W+,W- and center-of-mass energy
58 
59  IF(wp.LT.0.0 .OR. wm.LT.0.0) go to 1000
60 
61  IF(jout.EQ.0) THEN
62  IF(epp.LT.0.0) go to 1000
63  IF(epm.LT.0.0) go to 1000
64  IF(etp.LT.0.0) go to 1000
65  IF(etm.LT.0.0) go to 1000
66  IF(epp/(epm+0.01).LE.etp/(etm+0.01)) RETURN
67  ENDIF
68 C ********For strings which does not follow a jet-prod,
69 C scatter only if Ycm(JP)>Ycm(JT). When jets
70 C are produced just before this collision
71 C this requirement has already be enforced
72 C (see SUBROUTINE HIJHRD)
73  ihnt2(11)=jp
74  ihnt2(12)=jt
75 C
76 C
77 C
78  miss=0
79  pkc1=0.0
80  pkc2=0.0
81  pkc11=0.0
82  pkc12=0.0
83  pkc21=0.0
84  pkc22=0.0
85  dpkc11=0.0
86  dpkc12=0.0
87  dpkc21=0.0
88  dpkc22=0.0
89  IF(nfp(jp,10).EQ.1.OR.nft(jt,10).EQ.1) THEN
90  IF(nfp(jp,10).EQ.1) THEN
91  phi1=ulangl(pp(jp,10),pp(jp,11))
92  ppjet=sqrt(pp(jp,10)**2+pp(jp,11)**2)
93  pkc1=ppjet
94  pkc11=pp(jp,10)
95  pkc12=pp(jp,11)
96  ENDIF
97  IF(nft(jt,10).EQ.1) THEN
98  phi2=ulangl(pt(jt,10),pt(jt,11))
99  ptjet=sqrt(pt(jt,10)**2+pt(jt,11)**2)
100  pkc2=ptjet
101  pkc21=pt(jt,10)
102  pkc22=pt(jt,11)
103  ENDIF
104  IF(ihpr2(4).GT.0.AND.ihnt2(1).GT.1.AND.ihnt2(3).GT.1) THEN
105  IF(nfp(jp,10).EQ.0) THEN
106  phi=-phi2
107  ELSE IF(nft(jt,10).EQ.0) THEN
108  phi=phi1
109  ELSE
110  phi=(phi1+phi2-hipr1(40))/2.0
111  ENDIF
112  bx=hint1(19)*cos(hint1(20))
113  by=hint1(19)*sin(hint1(20))
114  xp0=yp(1,jp)
115  yp0=yp(2,jp)
116  xt0=yt(1,jt)+bx
117  yt0=yt(2,jt)+by
118  r1=max(1.2*ihnt2(1)**0.3333333,
119  & sqrt(xp0**2+yp0**2))
120  r2=max(1.2*ihnt2(3)**0.3333333,
121  & sqrt((xt0-bx)**2+(yt0-by)**2))
122  IF(abs(cos(phi)).LT.1.0e-5) THEN
123  dd1=r1
124  dd2=r1
125  dd3=abs(by+sqrt(r2**2-(xp0-bx)**2)-yp0)
126  dd4=abs(by-sqrt(r2**2-(xp0-bx)**2)-yp0)
127  go to 5
128  ENDIF
129  bb=2.0*sin(phi)*(cos(phi)*yp0-sin(phi)*xp0)
130  cc=(yp0**2-r1**2)*cos(phi)**2+xp0*sin(phi)*(
131  & xp0*sin(phi)-2.0*yp0*cos(phi))
132  dd=bb**2-4.0*cc
133  IF(dd.LT.0.0) go to 10
134  xx1=(-bb+sqrt(dd))/2.0
135  xx2=(-bb-sqrt(dd))/2.0
136  dd1=abs((xx1-xp0)/cos(phi))
137  dd2=abs((xx2-xp0)/cos(phi))
138 C
139  bb=2.0*sin(phi)*(cos(phi)*(yt0-by)-sin(phi)*xt0)-2.0*bx
140  cc=(bx**2+(yt0-by)**2-r2**2)*cos(phi)**2+xt0*sin(phi)
141  & *(xt0*sin(phi)-2.0*cos(phi)*(yt0-by))
142  & -2.0*bx*sin(phi)*(cos(phi)*(yt0-by)-sin(phi)*xt0)
143  dd=bb**2-4.0*cc
144  IF(dd.LT.0.0) go to 10
145  xx1=(-bb+sqrt(dd))/2.0
146  xx2=(-bb-sqrt(dd))/2.0
147  dd3=abs((xx1-xt0)/cos(phi))
148  dd4=abs((xx2-xt0)/cos(phi))
149 C
150  5 dd1=min(dd1,dd3)
151  dd2=min(dd2,dd4)
152  IF(dd1.LT.hipr1(13)) dd1=0.0
153  IF(dd2.LT.hipr1(13)) dd2=0.0
154  IF(nfp(jp,10).EQ.1.AND.ppjet.GT.hipr1(11)) THEN
155  dp1=dd1*hipr1(14)/2.0
156  dp1=min(dp1,ppjet-hipr1(11))
157  pkc1=ppjet-dp1
158  dpx1=cos(phi1)*dp1
159  dpy1=sin(phi1)*dp1
160  pkc11=pp(jp,10)-dpx1
161  pkc12=pp(jp,11)-dpy1
162  IF(dp1.GT.0.0) THEN
163  cthep=pp(jp,12)/sqrt(pp(jp,12)**2+ppjet**2)
164  dpz1=dp1*cthep/sqrt(1.0-cthep**2)
165  dpe1=sqrt(dpx1**2+dpy1**2+dpz1**2)
166  eppprm=pp(jp,4)+pp(jp,3)-dpe1-dpz1
167  epmprm=pp(jp,4)-pp(jp,3)-dpe1+dpz1
168  IF(eppprm.LE.0.0.OR.epmprm.LE.0.0) go to 15
169  epp=eppprm
170  epm=epmprm
171  pp(jp,10)=pkc11
172  pp(jp,11)=pkc12
173  npj(jp)=npj(jp)+1
174  kfpj(jp,npj(jp))=21
175  pjpx(jp,npj(jp))=dpx1
176  pjpy(jp,npj(jp))=dpy1
177  pjpz(jp,npj(jp))=dpz1
178  pjpe(jp,npj(jp))=dpe1
179  pjpm(jp,npj(jp))=0.0
180  pp(jp,3)=pp(jp,3)-dpz1
181  pp(jp,4)=pp(jp,4)-dpe1
182  ENDIF
183  ENDIF
184  15 IF(nft(jt,10).EQ.1.AND.ptjet.GT.hipr1(11)) THEN
185  dp2=dd2*hipr1(14)/2.0
186  dp2=min(dp2,ptjet-hipr1(11))
187  pkc2=ptjet-dp2
188  dpx2=cos(phi2)*dp2
189  dpy2=sin(phi2)*dp2
190  pkc21=pt(jt,10)-dpx2
191  pkc22=pt(jt,11)-dpy2
192  IF(dp2.GT.0.0) THEN
193  cthet=pt(jt,12)/sqrt(pt(jt,12)**2+ptjet**2)
194  dpz2=dp2*cthet/sqrt(1.0-cthet**2)
195  dpe2=sqrt(dpx2**2+dpy2**2+dpz2**2)
196  etpprm=pt(jt,4)+pt(jt,3)-dpe2-dpz2
197  etmprm=pt(jt,4)-pt(jt,3)-dpe2+dpz2
198  IF(etpprm.LE.0.0.OR.etmprm.LE.0.0) go to 16
199  etp=etpprm
200  etm=etmprm
201  pt(jt,10)=pkc21
202  pt(jt,11)=pkc22
203  ntj(jt)=ntj(jt)+1
204  kftj(jt,ntj(jt))=21
205  pjtx(jt,ntj(jt))=dpx2
206  pjty(jt,ntj(jt))=dpy2
207  pjtz(jt,ntj(jt))=dpz2
208  pjte(jt,ntj(jt))=dpe2
209  pjtm(jt,ntj(jt))=0.0
210  pt(jt,3)=pt(jt,3)-dpz2
211  pt(jt,4)=pt(jt,4)-dpe2
212  ENDIF
213  ENDIF
214  16 dpkc11=-(pp(jp,10)-pkc11)/2.0
215  dpkc12=-(pp(jp,11)-pkc12)/2.0
216  dpkc21=-(pt(jt,10)-pkc21)/2.0
217  dpkc22=-(pt(jt,11)-pkc22)/2.0
218  wp=epp+etp
219  wm=epm+etm
220  sw=wp*wm
221  ENDIF
222  ENDIF
223 C ********If jet is quenched the pt from valence quark
224 C hard scattering has to reduced by d*kapa
225 C
226 C
227 10 ptp02=pp(jp,1)**2+pp(jp,2)**2
228  ptt02=pt(jt,1)**2+pt(jt,2)**2
229 C
230  amq=max(pp(jp,14)+pp(jp,15),pt(jt,14)+pt(jt,15))
231  amx=hipr1(1)+amq
232 C ********consider mass cut-off for strings which
233 C must also include quark's mass
234  amp0=amx
235  dpm0=amx
236  nfdp=0
237  IF(nfp(jp,5).LE.2.AND.nfp(jp,3).NE.0) THEN
238  amp0=ulmass(nfp(jp,3))
239  nfdp=nfp(jp,3)+2*nfp(jp,3)/abs(nfp(jp,3))
240  dpm0=ulmass(nfdp)
241  IF(dpm0.LE.0.0) THEN
242  nfdp=nfdp-2*nfdp/abs(nfdp)
243  dpm0=ulmass(nfdp)
244  ENDIF
245  ENDIF
246  amt0=amx
247  dtm0=amx
248  nfdt=0
249  IF(nft(jt,5).LE.2.AND.nft(jt,3).NE.0) THEN
250  amt0=ulmass(nft(jt,3))
251  nfdt=nft(jt,3)+2*nft(jt,3)/abs(nft(jt,3))
252  dtm0=ulmass(nfdt)
253  IF(dtm0.LE.0.0) THEN
254  nfdt=nfdt-2*nfdt/abs(nfdt)
255  dtm0=ulmass(nfdt)
256  ENDIF
257  ENDIF
258 C
259  ampn=sqrt(amp0**2+ptp02)
260  amtn=sqrt(amt0**2+ptt02)
261  snn=(ampn+amtn)**2+0.001
262 C
263  IF(sw.LT.snn+0.001) go to 4000
264 C ********Scatter only if SW>SNN
265 C*****give some PT kick to the two exited strings******************
266 20 swptn=4.0*(max(amp0,amt0)**2+max(ptp02,ptt02))
267  swptd=4.0*(max(dpm0,dtm0)**2+max(ptp02,ptt02))
268  swptx=4.0*(amx**2+max(ptp02,ptt02))
269  IF(sw.LE.swptn) THEN
270  pkcmx=0.0
271  ELSE IF(sw.GT.swptn .AND. sw.LE.swptd
272  & .AND.npj(jp).EQ.0.AND.ntj(jt).EQ.0) THEN
273  pkcmx=sqrt(sw/4.0-max(amp0,amt0)**2)
274  & -sqrt(max(ptp02,ptt02))
275  ELSE IF(sw.GT.swptd .AND. sw.LE.swptx
276  & .AND.npj(jp).EQ.0.AND.ntj(jt).EQ.0) THEN
277  pkcmx=sqrt(sw/4.0-max(dpm0,dtm0)**2)
278  & -sqrt(max(ptp02,ptt02))
279  ELSE IF(sw.GT.swptx) THEN
280  pkcmx=sqrt(sw/4.0-amx**2)-sqrt(max(ptp02,ptt02))
281  ENDIF
282 C ********maximun PT kick
283 C*********************************************************
284 C
285  IF(nfp(jp,10).EQ.1.OR.nft(jt,10).EQ.1) THEN
286  IF(pkc1.GT.pkcmx) THEN
287  pkc1=pkcmx
288  pkc11=pkc1*cos(phi1)
289  pkc12=pkc1*sin(phi1)
290  dpkc11=-(pp(jp,10)-pkc11)/2.0
291  dpkc12=-(pp(jp,11)-pkc12)/2.0
292  ENDIF
293  IF(pkc2.GT.pkcmx) THEN
294  pkc2=pkcmx
295  pkc21=pkc2*cos(phi2)
296  pkc22=pkc2*sin(phi2)
297  dpkc21=-(pt(jt,10)-pkc21)/2.0
298  dpkc22=-(pt(jt,11)-pkc22)/2.0
299  ENDIF
300  dpkc1=dpkc11+dpkc21
301  dpkc2=dpkc12+dpkc22
302  nfp(jp,10)=-nfp(jp,10)
303  nft(jt,10)=-nft(jt,10)
304  go to 40
305  ENDIF
306 C ********If the valence quarks had a hard-collision
307 C the pt kick is the pt from hard-collision.
308  i_sng=0
309  IF(ihpr2(13).NE.0 .AND. atl_ran(nseed).LE.hidat(4)) i_sng=1
310  IF((nfp(jp,5).EQ.3 .OR.nft(jt,5).EQ.3).OR.
311  & (npj(jp).NE.0.OR.nfp(jp,10).NE.0).OR.
312  & (ntj(jt).NE.0.OR.nft(jt,10).NE.0)) i_sng=0
313 C
314 C ********decite whether to have single-diffractive
315  IF(ihpr2(5).EQ.0) THEN
316  pkc=hipr1(2)*sqrt(-alog(1.0-atl_ran(nseed)
317  & *(1.0-exp(-pkcmx**2/hipr1(2)**2))))
318  go to 30
319  ENDIF
320  pkc=hirnd2(3,0.0,pkcmx**2)
321  pkc=sqrt(pkc)
322  IF(pkc.GT.hipr1(20))
323  & pkc=hipr1(2)*sqrt(-alog(exp(-hipr1(20)**2/hipr1(2)**2)
324  & -atl_ran(nseed)*(exp(-hipr1(20)**2/hipr1(2)**2)-
325  & exp(-pkcmx**2/hipr1(2)**2))))
326 C
327  IF(i_sng.EQ.1) pkc=0.65*sqrt(
328  & -alog(1.0-atl_ran(nseed)*(1.0-exp(-pkcmx**2/0.65**2))))
329 C ********select PT kick
330 30 phi0=2.0*hipr1(40)*atl_ran(nseed)
331  pkc11=pkc*sin(phi0)
332  pkc12=pkc*cos(phi0)
333  pkc21=-pkc11
334  pkc22=-pkc12
335  dpkc1=0.0
336  dpkc2=0.0
337 40 pp11=pp(jp,1)+pkc11-dpkc1
338  pp12=pp(jp,2)+pkc12-dpkc2
339  pt11=pt(jt,1)+pkc21-dpkc1
340  pt12=pt(jt,2)+pkc22-dpkc2
341  ptp2=pp11**2+pp12**2
342  ptt2=pt11**2+pt12**2
343 C
344  ampn=sqrt(amp0**2+ptp2)
345  amtn=sqrt(amt0**2+ptt2)
346  snn=(ampn+amtn)**2+0.001
347 C***************************************
348  wp=epp+etp
349  wm=epm+etm
350  sw=wp*wm
351 C****************************************
352  IF(sw.LT.snn) THEN
353  miss=miss+1
354  IF(miss.LE.100) then
355  pkc=0.0
356  go to 30
357  ENDIF
358  IF(ihpr2(10).NE.0)
359  &WRITE(6,*) 'Error occured in Pt kick section of HIJSFT'
360  go to 4000
361  ENDIF
362 C******************************************************************
363  ampd=sqrt(dpm0**2+ptp2)
364  amtd=sqrt(dtm0**2+ptt2)
365 
366  ampx=sqrt(amx**2+ptp2)
367  amtx=sqrt(amx**2+ptt2)
368 
369  dpn=ampn**2/sw
370  dtn=amtn**2/sw
371  dpd=ampd**2/sw
372  dtd=amtd**2/sw
373  dpx=ampx**2/sw
374  dtx=amtx**2/sw
375 C
376  spntd=(ampn+amtd)**2
377  spntx=(ampn+amtx)**2
378 C ********CM energy if proj=N,targ=N*
379  spdtn=(ampd+amtn)**2
380  spxtn=(ampx+amtn)**2
381 C ********CM energy if proj=N*,targ=N
382  spdtx=(ampd+amtx)**2
383  spxtd=(ampx+amtd)**2
384  sdd=(ampd+amtd)**2
385  sxx=(ampx+amtx)**2
386 
387 C
388 C
389 C ********CM energy if proj=delta, targ=delta
390 C****************There are many different cases**********
391 c IF(IHPR2(15).EQ.1) GO TO 500
392 C
393 C ********to have DPM type soft interactions
394 C
395  45 CONTINUE
396  IF(sw.GT.sxx+0.001) THEN
397  IF(i_sng.EQ.0) THEN
398  d1=dpx
399  d2=dtx
400  nfp3=0
401  nft3=0
402  go to 400
403  ELSE
404 c**** 5/30/1998 this is identical to the above statement. Added to
405 c**** avoid questional branching to block.
406  IF((nfp(jp,5).EQ.3 .AND.nft(jt,5).EQ.3).OR.
407  & (npj(jp).NE.0.OR.nfp(jp,10).NE.0).OR.
408  & (ntj(jt).NE.0.OR.nft(jt,10).NE.0)) THEN
409  d1=dpx
410  d2=dtx
411  nfp3=0
412  nft3=0
413  go to 400
414  ENDIF
415 C ********do not allow excited strings to have
416 C single-diffr
417  IF(atl_ran(nseed).GT.0.5.OR.(nft(jt,5).GT.2.OR.
418  & ntj(jt).NE.0.OR.nft(jt,10).NE.0)) THEN
419  d1=dpn
420  d2=dtx
421  nfp3=nfp(jp,3)
422  nft3=0
423  go to 220
424  ELSE
425  d1=dpx
426  d2=dtn
427  nfp3=0
428  nft3=nft(jt,3)
429  go to 240
430  ENDIF
431 C ********have single diffractive collision
432  ENDIF
433  ELSE IF(sw.GT.max(spdtx,spxtd)+0.001 .AND.
434  & sw.LE.sxx+0.001) THEN
435  IF(((npj(jp).EQ.0.AND.ntj(jt).EQ.0.AND.
436  & atl_ran(nseed).GT.0.5).OR.(npj(jp).EQ.0
437  & .AND.ntj(jt).NE.0)).AND.nfp(jp,5).LE.2) THEN
438  d1=dpd
439  d2=dtx
440  nfp3=nfdp
441  nft3=0
442  go to 220
443  ELSE IF(ntj(jt).EQ.0.AND.nft(jt,5).LE.2) THEN
444  d1=dpx
445  d2=dtd
446  nfp3=0
447  nft3=nfdt
448  go to 240
449  ENDIF
450  go to 4000
451  ELSE IF(sw.GT.min(spdtx,spxtd)+0.001.AND.
452  & sw.LE.max(spdtx,spxtd)+0.001) THEN
453  IF(spdtx.LE.spxtd.AND.npj(jp).EQ.0
454  & .AND.nfp(jp,5).LE.2) THEN
455  d1=dpd
456  d2=dtx
457  nfp3=nfdp
458  nft3=0
459  go to 220
460  ELSE IF(spdtx.GT.spxtd.AND.ntj(jt).EQ.0
461  & .AND.nft(jt,5).LE.2) THEN
462  d1=dpx
463  d2=dtd
464  nfp3=0
465  nft3=nfdt
466  go to 240
467  ENDIF
468 c*** 5/30/1998 added to avoid questional branching to another block
469 c*** this is identical to the statement following the next ELSE IF
470  IF(((npj(jp).EQ.0.AND.ntj(jt).EQ.0
471  & .AND.atl_ran(nseed).GT.0.5).OR.(npj(jp).EQ.0
472  & .AND.ntj(jt).NE.0)).AND.nfp(jp,5).LE.2) THEN
473  d1=dpn
474  d2=dtx
475  nfp3=nfp(jp,3)
476  nft3=0
477  go to 220
478  ELSE IF(ntj(jt).EQ.0.AND.nft(jt,5).LE.2) THEN
479  d1=dpx
480  d2=dtn
481  nfp3=0
482  nft3=nft(jt,3)
483  go to 240
484  ENDIF
485  go to 4000
486  ELSE IF(sw.GT.max(spntx,spxtn)+0.001 .AND.
487  & sw.LE.min(spdtx,spxtd)+0.001) THEN
488  IF(((npj(jp).EQ.0.AND.ntj(jt).EQ.0
489  & .AND.atl_ran(nseed).GT.0.5).OR.(npj(jp).EQ.0
490  & .AND.ntj(jt).NE.0)).AND.nfp(jp,5).LE.2) THEN
491  d1=dpn
492  d2=dtx
493  nfp3=nfp(jp,3)
494  nft3=0
495  go to 220
496  ELSE IF(ntj(jt).EQ.0.AND.nft(jt,5).LE.2) THEN
497  d1=dpx
498  d2=dtn
499  nfp3=0
500  nft3=nft(jt,3)
501  go to 240
502  ENDIF
503  go to 4000
504  ELSE IF(sw.GT.min(spntx,spxtn)+0.001 .AND.
505  & sw.LE.max(spntx,spxtn)+0.001) THEN
506  IF(spntx.LE.spxtn.AND.npj(jp).EQ.0
507  & .AND.nfp(jp,5).LE.2) THEN
508  d1=dpn
509  d2=dtx
510  nfp3=nfp(jp,3)
511  nft3=0
512  go to 220
513  ELSEIF(spntx.GT.spxtn.AND.ntj(jt).EQ.0
514  & .AND.nft(jt,5).LE.2) THEN
515  d1=dpx
516  d2=dtn
517  nfp3=0
518  nft3=nft(jt,3)
519  go to 240
520  ENDIF
521  go to 4000
522  ELSE IF(sw.LE.min(spntx,spxtn)+0.001 .AND.
523  & (npj(jp).NE.0 .OR.ntj(jt).NE.0)) THEN
524  go to 4000
525  ELSE IF(sw.LE.min(spntx,spxtn)+0.001 .AND.
526  & nfp(jp,5).GT.2.AND.nft(jt,5).GT.2) THEN
527  go to 4000
528  ELSE IF(sw.GT.sdd+0.001.AND.sw.LE.
529  & min(spntx,spxtn)+0.001) THEN
530  d1=dpd
531  d2=dtd
532  nfp3=nfdp
533  nft3=nfdt
534  go to 100
535  ELSE IF(sw.GT.max(spntd,spdtn)+0.001
536  & .AND. sw.LE.sdd+0.001) THEN
537  IF(atl_ran(nseed).GT.0.5) THEN
538  d1=dpd
539  d2=dtn
540  nfp3=nfdp
541  nft3=nft(jt,3)
542  go to 100
543  ELSE
544  d1=dpn
545  d2=dtd
546  nfp3=nfp(jp,3)
547  nft3=nfdt
548  go to 100
549  ENDIF
550  ELSE IF(sw.GT.min(spntd,spdtn)+0.001
551  & .AND. sw.LE.max(spntd,spdtn)+0.001) THEN
552  IF(spntd.GT.spdtn) THEN
553  d1=dpd
554  d2=dtn
555  nfp3=nfdp
556  nft3=nft(jt,3)
557  go to 100
558  ELSE
559  d1=dpn
560  d2=dtd
561  nfp3=nfp(jp,3)
562  nft3=nfdt
563  go to 100
564  ENDIF
565  ELSE IF(sw.LE.min(spntd,spdtn)+0.001) THEN
566  d1=dpn
567  d2=dtn
568  nfp3=nfp(jp,3)
569  nft3=nft(jt,3)
570  go to 100
571  ENDIF
572  WRITE(6,*) 'Error in HIJSFT: There is no path to here'
573  RETURN
574 C
575 C*************** elastic scattering ***************
576 C this is like elastic, both proj and targ mass
577 C must be fixed
578 C***************************************************
579 100 nfp5=max(2,nfp(jp,5))
580  nft5=max(2,nft(jt,5))
581  bb1=1.0+d1-d2
582  bb2=1.0+d2-d1
583  IF(bb1**2.LT.4.0*d1 .OR. bb2**2.LT.4.0*d2) THEN
584  miss=miss+1
585  IF(miss.GT.100.OR.pkc.EQ.0.0) go to 3000
586  pkc=pkc*0.5
587  go to 30
588  ENDIF
589  IF(atl_ran(nseed).LT.0.5) THEN
590  x1=(bb1-sqrt(bb1**2-4.0*d1))/2.0
591  x2=(bb2-sqrt(bb2**2-4.0*d2))/2.0
592  ELSE
593  x1=(bb1+sqrt(bb1**2-4.0*d1))/2.0
594  x2=(bb2+sqrt(bb2**2-4.0*d2))/2.0
595  ENDIF
596  ihnt2(13)=2
597  go to 600
598 C
599 C********** Single diffractive ***********************
600 C either proj or targ's mass is fixed
601 C*****************************************************
602 220 nfp5=max(2,nfp(jp,5))
603  nft5=3
604  IF(nfp3.EQ.0) nfp5=3
605  bb2=1.0+d2-d1
606  IF(bb2**2.LT.4.0*d2) THEN
607  miss=miss+1
608  IF(miss.GT.100.OR.pkc.EQ.0.0) go to 3000
609  pkc=pkc*0.5
610  go to 30
611  ENDIF
612  xmin=(bb2-sqrt(bb2**2-4.0*d2))/2.0
613  xmax=(bb2+sqrt(bb2**2-4.0*d2))/2.0
614  miss4=0
615 222 x2=hirnd2(6,xmin,xmax)
616  x1=d1/(1.0-x2)
617  IF(x2*(1.0-x1).LT.(d2+1.e-4/sw)) THEN
618  miss4=miss4+1
619  IF(miss4.LE.1000) go to 222
620  go to 5000
621  ENDIF
622  ihnt2(13)=2
623  go to 600
624 C ********Fix proj mass*********
625 240 nfp5=3
626  nft5=max(2,nft(jt,5))
627  IF(nft3.EQ.0) nft5=3
628  bb1=1.0+d1-d2
629  IF(bb1**2.LT.4.0*d1) THEN
630  miss=miss+1
631  IF(miss.GT.100.OR.pkc.EQ.0.0) go to 3000
632  pkc=pkc*0.5
633  go to 30
634  ENDIF
635  xmin=(bb1-sqrt(bb1**2-4.0*d1))/2.0
636  xmax=(bb1+sqrt(bb1**2-4.0*d1))/2.0
637  miss4=0
638 242 x1=hirnd2(6,xmin,xmax)
639  x2=d2/(1.0-x1)
640  IF(x1*(1.0-x2).LT.(d1+1.e-4/sw)) THEN
641  miss4=miss4+1
642  IF(miss4.LE.1000) go to 242
643  go to 5000
644  ENDIF
645  ihnt2(13)=2
646  go to 600
647 C ********Fix targ mass*********
648 C
649 C*************non-single diffractive**********************
650 C both proj and targ may not be fixed in mass
651 C*********************************************************
652 C
653 400 nfp5=3
654  nft5=3
655  bb1=1.0+d1-d2
656  bb2=1.0+d2-d1
657  IF(bb1**2.LT.4.0*d1 .OR. bb2**2.LT.4.0*d2) THEN
658  miss=miss+1
659  IF(miss.GT.100.OR.pkc.EQ.0.0) go to 3000
660  pkc=pkc*0.5
661  go to 30
662  ENDIF
663  xmin1=(bb1-sqrt(bb1**2-4.0*d1))/2.0
664  xmax1=(bb1+sqrt(bb1**2-4.0*d1))/2.0
665  xmin2=(bb2-sqrt(bb2**2-4.0*d2))/2.0
666  xmax2=(bb2+sqrt(bb2**2-4.0*d2))/2.0
667  miss4=0
668 410 x1=hirnd2(4,xmin1,xmax1)
669  x2=hirnd2(4,xmin2,xmax2)
670  IF(nfp(jp,5).EQ.3.OR.nft(jt,5).EQ.3) THEN
671  x1=hirnd2(6,xmin1,xmax1)
672  x2=hirnd2(6,xmin2,xmax2)
673  ENDIF
674 C ********
675  IF(abs(nfp(jp,1)*nfp(jp,2)).GT.1000000.OR.
676  & abs(nfp(jp,1)*nfp(jp,2)).LT.100) THEN
677  x1=hirnd2(5,xmin1,xmax1)
678  ENDIF
679  IF(abs(nft(jt,1)*nft(jt,2)).GT.1000000.OR.
680  & abs(nft(jt,1)*nft(jt,2)).LT.100) THEN
681  x2=hirnd2(5,xmin2,xmax2)
682  ENDIF
683 c IF(IOPMAIN.EQ.3) X1=HIRND2(6,XMIN1,XMAX1)
684 c IF(IOPMAIN.EQ.2) X2=HIRND2(6,XMIN2,XMAX2)
685 C ********For q-qbar or (qq)-(qq)bar system use symetric
686 C distribution, for q-(qq) or qbar-(qq)bar use
687 C unsymetrical distribution
688 C
689  IF(abs(nfp(jp,1)*nfp(jp,2)).GT.1000000) x1=1.0-x1
690  xxp=x1*(1.0-x2)
691  xxt=x2*(1.0-x1)
692  IF(xxp.LT.(d1+1.e-4/sw) .OR. xxt.LT.(d2+1.e-4/sw)) THEN
693  miss4=miss4+1
694  IF(miss4.LE.1000) go to 410
695  go to 5000
696  ENDIF
697  ihnt2(13)=3
698 C***************************************************
699 C***************************************************
700 600 CONTINUE
701  IF(x1*(1.0-x2).LT.(ampn**2-1.e-4)/sw.OR.
702  & x2*(1.0-x1).LT.(amtn**2-1.e-4)/sw) THEN
703  miss=miss+1
704  IF(miss.GT.100.OR.pkc.EQ.0.0) go to 2000
705  pkc=0.0
706  go to 30
707  ENDIF
708 C
709  epp=(1.0-x2)*wp
710  epm=x1*wm
711  etp=x2*wp
712  etm=(1.0-x1)*wm
713  pp(jp,3)=(epp-epm)/2.0
714  pp(jp,4)=(epp+epm)/2.0
715  IF(epp*epm-ptp2.LT.0.0) go to 6000
716  pp(jp,5)=sqrt(epp*epm-ptp2)
717  nfp(jp,3)=nfp3
718  nfp(jp,5)=nfp5
719 
720  pt(jt,3)=(etp-etm)/2.0
721  pt(jt,4)=(etp+etm)/2.0
722  IF(etp*etm-ptt2.LT.0.0) go to 6000
723  pt(jt,5)=sqrt(etp*etm-ptt2)
724  nft(jt,3)=nft3
725  nft(jt,5)=nft5
726 C*****recoil PT from hard-inter is shared by two end-partons
727 C so that pt=p1+p2
728  pp(jp,1)=pp11-pkc11
729  pp(jp,2)=pp12-pkc12
730 
731  kickdip=1
732  kickdit=1
733  IF(abs(nfp(jp,1)*nfp(jp,2)).GT.1000000.OR.
734  & abs(nfp(jp,1)*nfp(jp,2)).LT.100) THEN
735  kickdip=0
736  ENDIF
737  IF(abs(nft(jt,1)*nft(jt,2)).GT.1000000.OR.
738  & abs(nft(jt,1)*nft(jt,2)).LT.100) THEN
739  kickdit=0
740  ENDIF
741  IF((kickdip.EQ.0.AND.atl_ran(nseed).LT.0.5)
742  & .OR.(kickdip.NE.0.AND.atl_ran(nseed)
743  & .LT.0.5/(1.0+(pkc11**2+pkc12**2)/hipr1(22)**2))) THEN
744  pp(jp,6)=(pp(jp,1)-pp(jp,6)-pp(jp,8)-dpkc1)/2.0+pp(jp,6)
745  pp(jp,7)=(pp(jp,2)-pp(jp,7)-pp(jp,9)-dpkc2)/2.0+pp(jp,7)
746  pp(jp,8)=(pp(jp,1)-pp(jp,6)-pp(jp,8)-dpkc1)/2.0
747  & +pp(jp,8)+pkc11
748  pp(jp,9)=(pp(jp,2)-pp(jp,7)-pp(jp,9)-dpkc2)/2.0
749  & +pp(jp,9)+pkc12
750  ELSE
751  pp(jp,8)=(pp(jp,1)-pp(jp,6)-pp(jp,8)-dpkc1)/2.0+pp(jp,8)
752  pp(jp,9)=(pp(jp,2)-pp(jp,7)-pp(jp,9)-dpkc2)/2.0+pp(jp,9)
753  pp(jp,6)=(pp(jp,1)-pp(jp,6)-pp(jp,8)-dpkc1)/2.0
754  & +pp(jp,6)+pkc11
755  pp(jp,7)=(pp(jp,2)-pp(jp,7)-pp(jp,9)-dpkc2)/2.0
756  & +pp(jp,7)+pkc12
757  ENDIF
758  pp(jp,1)=pp(jp,6)+pp(jp,8)
759  pp(jp,2)=pp(jp,7)+pp(jp,9)
760 C ********pt kick for proj
761  pt(jt,1)=pt11-pkc21
762  pt(jt,2)=pt12-pkc22
763  IF((kickdit.EQ.0.AND.atl_ran(nseed).LT.0.5)
764  & .OR.(kickdit.NE.0.AND.atl_ran(nseed)
765  & .LT.0.5/(1.0+(pkc21**2+pkc22**2)/hipr1(22)**2))) THEN
766  pt(jt,6)=(pt(jt,1)-pt(jt,6)-pt(jt,8)-dpkc1)/2.0+pt(jt,6)
767  pt(jt,7)=(pt(jt,2)-pt(jt,7)-pt(jt,9)-dpkc2)/2.0+pt(jt,7)
768  pt(jt,8)=(pt(jt,1)-pt(jt,6)-pt(jt,8)-dpkc1)/2.0
769  & +pt(jt,8)+pkc21
770  pt(jt,9)=(pt(jt,2)-pt(jt,7)-pt(jt,9)-dpkc2)/2.0
771  & +pt(jt,9)+pkc22
772  ELSE
773  pt(jt,8)=(pt(jt,1)-pt(jt,6)-pt(jt,8)-dpkc1)/2.0+pt(jt,8)
774  pt(jt,9)=(pt(jt,2)-pt(jt,7)-pt(jt,9)-dpkc2)/2.0+pt(jt,9)
775  pt(jt,6)=(pt(jt,1)-pt(jt,6)-pt(jt,8)-dpkc1)/2.0
776  & +pt(jt,6)+pkc21
777  pt(jt,7)=(pt(jt,2)-pt(jt,7)-pt(jt,9)-dpkc2)/2.0
778  & +pt(jt,7)+pkc22
779  ENDIF
780  pt(jt,1)=pt(jt,6)+pt(jt,8)
781  pt(jt,2)=pt(jt,7)+pt(jt,9)
782 C ********pt kick for targ
783 
784  IF(npj(jp).NE.0) nfp(jp,5)=3
785  IF(ntj(jt).NE.0) nft(jt,5)=3
786 C ********jets must be connected to string
787  IF(epp/(epm+0.0001).LT.etp/(etm+0.0001).AND.
788  & abs(nfp(jp,1)*nfp(jp,2)).LT.1000000)THEN
789  DO 620 jsb=1,15
790  psb=pp(jp,jsb)
791  pp(jp,jsb)=pt(jt,jsb)
792  pt(jt,jsb)=psb
793  nsb=nfp(jp,jsb)
794  nfp(jp,jsb)=nft(jt,jsb)
795  nft(jt,jsb)=nsb
796 620 CONTINUE
797 C ********when Ycm(JP)<Ycm(JT) after the collision
798 C exchange the positions of the two
799  ENDIF
800 C
801  RETURN
802 C**************************************************
803 C**************************************************
804 1000 ierror=1
805  IF(ihpr2(10).EQ.0) RETURN
806  WRITE(6,*) ' Fatal HIJSFT start error,abandon this event'
807  WRITE(6,*) ' PROJ E+,E-,W+',epp,epm,wp
808  WRITE(6,*) ' TARG E+,E-,W-',etp,etm,wm
809  WRITE(6,*) ' W+*W-, (APN+ATN)^2',sw,snn
810  RETURN
811 2000 ierror=0
812  IF(ihpr2(10).EQ.0) RETURN
813  WRITE(6,*) ' (2)energy partition fail,'
814  WRITE(6,*) ' HIJSFT not performed, but continue'
815  WRITE(6,*) ' MP1,MPN',x1*(1.0-x2)*sw,ampn**2
816  WRITE(6,*) ' MT2,MTN',x2*(1.0-x1)*sw,amtn**2
817  RETURN
818 3000 ierror=0
819  IF(ihpr2(10).EQ.0) RETURN
820  WRITE(6,*) ' (3)something is wrong with the pt kick, '
821  WRITE(6,*) ' HIJSFT not performed, but continue'
822  WRITE(6,*) ' D1=',d1,' D2=',d2,' SW=',sw
823  WRITE(6,*) ' HISTORY NFP5=',nfp(jp,5),' NFT5=',nft(jt,5)
824  WRITE(6,*) ' THIS COLLISON NFP5=',nfp5, ' NFT5=',nft5
825  WRITE(6,*) ' # OF JET IN PROJ',npj(jp),' IN TARG',ntj(jt)
826  RETURN
827 4000 ierror=0
828  IF(ihpr2(10).EQ.0) RETURN
829  WRITE(6,*) ' (4)unable to choose process, but not harmful'
830  WRITE(6,*) ' HIJSFT not performed, but continue'
831  WRITE(6,*) ' PTP=',sqrt(ptp2),' PTT=',sqrt(ptt2),' SW=',sw
832  WRITE(6,*) ' AMCUT=',amx,' JP=',jp,' JT=',jt
833  WRITE(6,*) ' HISTORY NFP5=',nfp(jp,5),' NFT5=',nft(jt,5)
834  RETURN
835 5000 ierror=0
836  IF(ihpr2(10).EQ.0) RETURN
837  WRITE(6,*) ' energy partition failed(5),for limited try'
838  WRITE(6,*) ' HIJSFT not performed, but continue'
839  WRITE(6,*) ' NFP5=',nfp5,' NFT5=',nft5
840  WRITE(6,*) ' D1',d1,' X1(1-X2)',x1*(1.0-x2)
841  WRITE(6,*) ' D2',d2,' X2(1-X1)',x2*(1.0-x1)
842  RETURN
843 6000 pkc=0.0
844  miss=miss+1
845  IF(miss.LT.100) go to 30
846  ierror=1
847  IF(ihpr2(10).EQ.0) RETURN
848  WRITE(6,*) 'ERROR OCCURED, HIJSFT NOT PERFORMED'
849  WRITE(6,*) ' Abort this event'
850  WRITE(6,*) 'MTP,PTP2',epp*epm,ptp2,' MTT,PTT2',etp*etm,ptt2
851  RETURN
852  END