3 C*********************************************************************
5 SUBROUTINE PYSSPA_HIJING(IPU1,IPU2)
7 C...Generates spacelike parton showers.
8 IMPLICIT DOUBLE PRECISION(D)
9 #include "lujets_hijing.inc"
10 #include "ludat1_hijing.inc"
11 #include "ludat2_hijing.inc"
12 #include "pysubs_hijing.inc"
13 #include "pypars_hijing.inc"
14 #include "pyint1_hijing.inc"
15 #include "pyint2_hijing.inc"
16 #include "pyint3_hijing.inc"
17 DIMENSION KFLS(4),IS(2),XS(2),ZS(2),Q2S(2),TEVS(2),ROBO(5),
18 &XFS(2,-6:6),XFA(-6:6),XFB(-6:6),XFN(-6:6),WTAP(-6:6),WTSF(-6:6),
19 &THE2(2),ALAM(2),DQ2(3),DPC(3),DPD(4),DPB(4)
21 C...Calculate maximum virtuality and check that evolution possible.
26 IF(ISET(ISUB).EQ.1) THEN
28 ELSEIF(ISET(ISUB).EQ.3.OR.ISET(ISUB).EQ.4) THEN
30 IF(ISUB.EQ.8.OR.ISUB.EQ.76.OR.ISUB.EQ.77) Q2E=PMAS(24,1)**2
32 TMAX=LOG(PARP(67)*PARP(63)*Q2E/PARP(61)**2)
33 IF(PARP(67)*Q2E.LT.MAX(PARP(62)**2,2.*PARP(61)**2).OR.
36 C...Common constants and initial values. Save normal Lambda value.
37 XE0=2.*PARP(65)/VINT(1)
52 110 XFS(JT,KFL)=XSFX(JT,KFL)
54 IF(ISET(ISUB).EQ.3.OR.ISET(ISUB).EQ.4) DSH=VINT(26)*VINT(2)
56 C...Pick up leg with highest virtuality.
59 IF(N.GT.NS+1.AND.Q2S(2).GT.Q2S(1)) JT=2
63 130 XFB(KFL)=XFS(JT,KFL)
66 XE=MAX(XE0,XB*(1./(1.-PARP(66))-1.))
67 IF(XB+XE.GE.0.999) THEN
72 C...Maximum Q2 without or with Q2 ordering. Effective Lambda and n_f.
73 IF(MSTP(62).LE.1) THEN
74 Q2B=0.5*(1./ZS(JT)+1.)*Q2S(JT)+0.5*(1./ZS(JT)-1.)*(Q2S(3-JT)-
75 & SNGL(DSH)+SQRT((SNGL(DSH)+Q2S(1)+Q2S(2))**2+8.*Q2S(1)*Q2S(2)*
76 & ZS(JT)/(1.-ZS(JT))))
77 TEVB=LOG(PARP(63)*Q2B/ALAM(JT)**2)
82 ALSDUM=ULALPS_HIJING(PARP(63)*Q2B)
83 TEVB=TEVB+2.*LOG(ALAM(JT)/PARU(117))
86 B0=(33.-2.*MSTU(118))/6.
88 C...Calculate Altarelli-Parisi and structure function weights.
93 WTAPQ=16.*(1.-SQRT(XB+XE))/(3.*SQRT(XB))
94 DO 150 KFL=-MSTP(54),MSTP(54)
95 IF(KFL.EQ.0) WTAP(KFL)=6.*LOG((1.-XB)/XE)
96 150 IF(KFL.NE.0) WTAP(KFL)=WTAPQ
98 WTAP(0)=0.5*XB*(1./(XB+XE)-1.)
99 WTAP(KFLB)=8.*LOG((1.-XB)*(XB+XE)/XE)/3.
102 IF(KFLB.NE.21) XFBO=XFB(KFLB)
103 IF(KFLB.EQ.21) XFBO=XFB(0)
104 C***************************************************************
105 C**********ERROR HAS OCCURED HERE
108 WRITE(MSTU(11),1001) KFLB,XFB(KFLB)
111 C****************************************************************
112 DO 170 KFL=-MSTP(54),MSTP(54)
113 WTSF(KFL)=XFB(KFL)/XFBO
114 170 WTSUM=WTSUM+WTAP(KFL)*WTSF(KFL)
115 WTSUM=MAX(0.0001,WTSUM)
117 C...Choose new t: fix alpha_s, alpha_s(Q2), alpha_s(k_T2).
118 180 IF(MSTP(64).LE.0) THEN
119 TEVB=TEVB+LOG(RLU_HIJING(0))*PARU(2)/(PARU(111)*WTSUM)
120 ELSEIF(MSTP(64).EQ.1) THEN
121 TEVB=TEVB*EXP(MAX(-100.,LOG(RLU_HIJING(0))*B0/WTSUM))
123 TEVB=TEVB*EXP(MAX(-100.,LOG(RLU_HIJING(0))*B0/(5.*WTSUM)))
125 190 Q2REF=ALAM(JT)**2*EXP(TEVB)
128 C...Evolution ended or select flavour for branching parton.
129 IF(Q2B.LT.PARP(62)**2) THEN
132 WTRAN=RLU_HIJING(0)*WTSUM
135 WTRAN=WTRAN-WTAP(KFLA)*WTSF(KFLA)
136 IF(KFLA.LT.MSTP(54).AND.WTRAN.GT.0.) GOTO 200
137 IF(KFLA.EQ.0) KFLA=21
139 C...Choose z value and corrective weight.
140 IF(KFLB.EQ.21.AND.KFLA.EQ.21) THEN
141 Z=1./(1.+((1.-XB)/XB)*(XE/(1.-XB))**RLU_HIJING(0))
143 ELSEIF(KFLB.EQ.21) THEN
144 Z=XB/(1.-RLU_HIJING(0)*(1.-SQRT(XB+XE)))**2
145 WTZ=0.5*(1.+(1.-Z)**2)*SQRT(Z)
146 ELSEIF(KFLA.EQ.21) THEN
147 Z=XB*(1.+RLU_HIJING(0)*(1./(XB+XE)-1.))
150 Z=1.-(1.-XB)*(XE/((XB+XE)*(1.-XB)))**RLU_HIJING(0)
154 C...Option with resummation of soft gluon emission as effective z shift.
155 IF(MSTP(65).GE.1) THEN
157 IF(KFLB.NE.21) RSOFT=8./3.
158 Z=Z*(TEVB/TEVS(JT))**(RSOFT*XE/((XB+XE)*B0))
162 C...Option with alpha_s(k_T2)Q2): demand k_T2 > cutoff, reweight.
163 IF(MSTP(64).GE.2) THEN
164 IF((1.-Z)*Q2B.LT.PARP(62)**2) GOTO 180
165 ALPRAT=TEVB/(TEVB+LOG(1.-Z))
166 IF(ALPRAT.LT.5.*RLU_HIJING(0)) GOTO 180
167 IF(ALPRAT.GT.5.) WTZ=WTZ*ALPRAT/5.
170 C...Option with angular ordering requirement.
171 IF(MSTP(62).GE.3) THEN
172 THE2T=(4.*Z**2*Q2B)/(VINT(2)*(1.-Z)*XB**2)
173 IF(THE2T.GT.THE2(JT)) GOTO 180
176 C...Weighting with new structure functions.
177 CALL PYSTFU_HIJING(MINT(10+JT),XB,Q2REF,XFN,JT)
178 IF(KFLB.NE.21) XFBN=XFN(KFLB)
179 IF(KFLB.EQ.21) XFBN=XFN(0)
180 IF(XFBN.LT.1E-20) THEN
181 IF(KFLA.EQ.KFLB) THEN
185 ELSEIF(TEVBSV-TEVB.GT.0.2) THEN
186 TEVB=0.5*(TEVBSV+TEVB)
192 DO 210 KFL=-MSTP(54),MSTP(54)
193 210 XFB(KFL)=XFN(KFL)
195 CALL PYSTFU_HIJING(MINT(10+JT),XA,Q2REF,XFA,JT)
196 IF(KFLA.NE.21) XFAN=XFA(KFLA)
197 IF(KFLA.EQ.21) XFAN=XFA(0)
198 IF(XFAN.LT.1E-20) GOTO 160
199 IF(KFLA.NE.21) WTSFA=WTSF(KFLA)
200 IF(KFLA.EQ.21) WTSFA=WTSF(0)
201 IF(WTZ*XFAN/XFBN.LT.RLU_HIJING(0)*WTSFA) GOTO 160
204 C...Define two hard scatterers in their CM-frame.
205 220 IF(N.EQ.NS+2) THEN
207 DPLCM=SQRT((DSH+DQ2(1)+DQ2(2))**2-4D0*DQ2(1)*DQ2(2))/DSHR
210 IF(JR.EQ.1) IPO=IPUS1
211 IF(JR.EQ.2) IPO=IPUS2
220 P(I,3)=DPLCM*(-1)**(JR+1)
221 P(I,4)=(DSH+DQ2(3-JR)-DQ2(JR))/DSHR
222 P(I,5)=-SQRT(SNGL(DQ2(JR)))
225 K(IPO,4)=MOD(K(IPO,4),MSTU(5))+MSTU(5)*I
226 240 K(IPO,5)=MOD(K(IPO,5),MSTU(5))+MSTU(5)*I
228 C...Find maximum allowed mass of timelike parton.
229 ELSEIF(N.GT.NS+2) THEN
234 DPC(3)=0.5*(ABS(P(IS(1),3))+ABS(P(IS(2),3)))
235 DPD(1)=DSH+DQ2(JR)+DQ2(JT)
236 DPD(2)=DSHZ+DQ2(JR)+DQ2(3)
237 DPD(3)=SQRT(DPD(1)**2-4D0*DQ2(JR)*DQ2(JT))
238 DPD(4)=SQRT(DPD(2)**2-4D0*DQ2(JR)*DQ2(3))
240 IF(Q2S(JR).GE.(0.5*PARP(62))**2.AND.DPD(1)-DPD(3).GE.
241 & 1D-10*DPD(1)) IKIN=1
242 IF(IKIN.EQ.0) DMSMA=(DQ2(JT)/DBLE(ZS(JT))-DQ2(3))*(DSH/
243 & (DSH+DQ2(JT))-DSH/(DSHZ+DQ2(3)))
244 IF(IKIN.EQ.1) DMSMA=(DPD(1)*DPD(2)-DPD(3)*DPD(4))/(2.*
245 & DQ2(JR))-DQ2(JT)-DQ2(3)
247 C...Generate timelike parton shower (if required).
255 IF(KFLB.EQ.21.AND.KFLS(JT+2).NE.21) K(IT,2)=-KFLS(JT+2)
256 IF(KFLB.NE.21.AND.KFLS(JT+2).EQ.21) K(IT,2)=KFLB
257 P(IT,5)=ULMASS_HIJING(K(IT,2))
258 IF(SNGL(DMSMA).LE.P(IT,5)**2) GOTO 100
259 IF(MSTP(63).GE.1) THEN
260 P(IT,4)=(DSHZ-DSH-P(IT,5)**2)/DSHR
261 P(IT,3)=SQRT(P(IT,4)**2-P(IT,5)**2)
262 IF(MSTP(63).EQ.1) THEN
264 ELSEIF(MSTP(63).EQ.2) THEN
265 Q2TIM=MIN(SNGL(DMSMA),PARP(71)*Q2S(JT))
267 C'''Here remains to introduce angular ordering in first branching.
270 CALL LUSHOW_HIJING(IT,0,SQRT(Q2TIM))
271 IF(N.GE.IT+1) P(IT,5)=P(IT+1,5)
274 C...Reconstruct kinematics of branching: timelike parton shower.
276 IF(IKIN.EQ.0) DPT2=(DMSMA-DMS)*(DSHZ+DQ2(3))/(DSH+DQ2(JT))
277 IF(IKIN.EQ.1) DPT2=(DMSMA-DMS)*(0.5*DPD(1)*DPD(2)+0.5*DPD(3)*
278 & DPD(4)-DQ2(JR)*(DQ2(JT)+DQ2(3)+DMS))/(4.*DSH*DPC(3)**2)
279 IF(DPT2.LT.0.) GOTO 100
280 DPB(1)=(0.5*DPD(2)-DPC(JR)*(DSHZ+DQ2(JR)-DQ2(JT)-DMS)/
281 & DSHR)/DPC(3)-DPC(3)
282 P(IT,1)=SQRT(SNGL(DPT2))
283 P(IT,3)=DPB(1)*(-1)**(JT+1)
284 P(IT,4)=(DSHZ-DSH-DMS)/DSHR
286 DPB(1)=SQRT(DPB(1)**2+DPT2)
287 DPB(2)=SQRT(DPB(1)**2+DMS)
289 DPB(4)=SQRT(DPB(3)**2+DMS)
290 DBEZ=(DPB(4)*DPB(1)-DPB(3)*DPB(2))/(DPB(4)*DPB(2)-DPB(3)*
292 CALL LUDBRB_HIJING(IT+1,N,0.,0.,0D0,0D0,DBEZ)
293 THE=ULANGL_HIJING(P(IT,3),P(IT,1))
294 CALL LUDBRB_HIJING(IT+1,N,THE,0.,0D0,0D0,0D0)
297 C...Reconstruct kinematics of branching: spacelike parton.
305 P(N+1,3)=P(IT,3)+P(IS(JT),3)
306 P(N+1,4)=P(IT,4)+P(IS(JT),4)
307 P(N+1,5)=-SQRT(SNGL(DQ2(3)))
309 C...Define colour flow of branching.
313 IF((K(N+1,2).GT.0.AND.K(N+1,2).NE.21.AND.K(ID1,2).GT.0.AND.
314 & K(ID1,2).NE.21).OR.(K(N+1,2).LT.0.AND.K(ID1,2).EQ.21).OR.
315 & (K(N+1,2).EQ.21.AND.K(ID1,2).EQ.21.AND.RLU_HIJING(0).GT.0.5).OR.
316 & (K(N+1,2).EQ.21.AND.K(ID1,2).LT.0)) ID1=IS(JT)
318 K(N+1,4)=K(N+1,4)+ID1
319 K(N+1,5)=K(N+1,5)+ID2
320 K(ID1,4)=K(ID1,4)+MSTU(5)*(N+1)
321 K(ID1,5)=K(ID1,5)+MSTU(5)*ID2
322 K(ID2,4)=K(ID2,4)+MSTU(5)*ID1
323 K(ID2,5)=K(ID2,5)+MSTU(5)*(N+1)
326 C...Boost to new CM-frame.
327 CALL LUDBRB_HIJING(NS+1,N,0.,0.,-DBLE((P(N,1)+P(IS(JR),1))/(P(N
328 $ ,4)+P(IS(JR),4))),0D0,-DBLE((P(N,3)+P(IS(JR),3))/(P(N,4)
330 IR=N+(JT-1)*(IS(1)-N)
331 CALL LUDBRB_HIJING(NS+1,N,-ULANGL_HIJING(P(IR,3),P(IR,1)),PARU(2
332 $ )*RLU_HIJING(0),0D0,0D0,0D0)
335 C...Save quantities, loop back.
339 IF(MSTP(62).GE.3) THE2(JT)=THE2T
341 IF(Q2B.GE.(0.5*PARP(62))**2) THEN
347 270 XFS(JT,KFL)=XFA(KFL)
353 IF(N.GT.MSTU(4)-MSTU(32)-10) THEN
354 CALL LUERRM_HIJING(11
355 $ ,'(PYSSPA_HIJING:) no more memory left in LUJETS_HIJING')
356 IF(MSTU(21).GE.1) N=NS
357 IF(MSTU(21).GE.1) RETURN
359 IF(MAX(Q2S(1),Q2S(2)).GE.(0.5*PARP(62))**2.OR.N.LE.NS+1) GOTO 120
361 C...Boost hard scattering partons to frame of shower initiators.
363 280 ROBO(J+2)=(P(NS+1,J)+P(NS+2,J))/(P(NS+1,4)+P(NS+2,4))
365 290 P(N+2,J)=P(NS+1,J)
366 ROBOT=ROBO(3)**2+ROBO(4)**2+ROBO(5)**2
367 IF(ROBOT.GE.0.999999) THEN
368 ROBOT=1.00001*SQRT(ROBOT)
369 ROBO(3)=ROBO(3)/ROBOT
370 ROBO(4)=ROBO(4)/ROBOT
371 ROBO(5)=ROBO(5)/ROBOT
373 CALL LUDBRB_HIJING(N+2,N+2,0.,0.,-DBLE(ROBO(3)),-DBLE(ROBO(4)),
375 ROBO(2)=ULANGL_HIJING(P(N+2,1),P(N+2,2))
376 ROBO(1)=ULANGL_HIJING(P(N+2,3),SQRT(P(N+2,1)**2+P(N+2,2)**2))
377 CALL LUDBRB_HIJING(MINT(83)+5,NS,ROBO(1),ROBO(2),DBLE(ROBO(3)),
378 &DBLE(ROBO(4)),DBLE(ROBO(5)))
380 C...Store user information. Reset Lambda value.
385 300 VINT(140+JT)=XS(JT)
387 1000 FORMAT(5X,'structure function has a zero point here')
388 1001 FORMAT(5X,'xf(x,i=',I5,')=',F10.5)