]> git.uio.no Git - u/mrichter/AliRoot.git/blame - HIJING/hipyset1_35/pysspa_hijing.F
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
fWSN->Eval(0.001) to avoid fpe.
[u/mrichter/AliRoot.git] / HIJING / hipyset1_35 / pysspa_hijing.F
CommitLineData
e74335a4 1* $Id$
2
3C*********************************************************************
4
5 SUBROUTINE PYSSPA_HIJING(IPU1,IPU2)
6
7C...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)
20
21C...Calculate maximum virtuality and check that evolution possible.
22 IPUS1=IPU1
23 IPUS2=IPU2
24 ISUB=MINT(1)
25 Q2E=VINT(52)
26 IF(ISET(ISUB).EQ.1) THEN
27 Q2E=Q2E/PARP(67)
28 ELSEIF(ISET(ISUB).EQ.3.OR.ISET(ISUB).EQ.4) THEN
29 Q2E=PMAS(23,1)**2
30 IF(ISUB.EQ.8.OR.ISUB.EQ.76.OR.ISUB.EQ.77) Q2E=PMAS(24,1)**2
31 ENDIF
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.
34 &TMAX.LT.0.2) RETURN
35
36C...Common constants and initial values. Save normal Lambda value.
37 XE0=2.*PARP(65)/VINT(1)
38 ALAMS=PARU(111)
39 PARU(111)=PARP(61)
40 NS=N
41 100 N=NS
42 DO 110 JT=1,2
43 KFLS(JT)=MINT(14+JT)
44 KFLS(JT+2)=KFLS(JT)
45 XS(JT)=VINT(40+JT)
46 ZS(JT)=1.
47 Q2S(JT)=PARP(67)*Q2E
48 TEVS(JT)=TMAX
49 ALAM(JT)=PARP(61)
50 THE2(JT)=100.
51 DO 110 KFL=-6,6
52 110 XFS(JT,KFL)=XSFX(JT,KFL)
53 DSH=VINT(44)
54 IF(ISET(ISUB).EQ.3.OR.ISET(ISUB).EQ.4) DSH=VINT(26)*VINT(2)
55
56C...Pick up leg with highest virtuality.
57 120 N=N+1
58 JT=1
59 IF(N.GT.NS+1.AND.Q2S(2).GT.Q2S(1)) JT=2
60 KFLB=KFLS(JT)
61 XB=XS(JT)
62 DO 130 KFL=-6,6
63 130 XFB(KFL)=XFS(JT,KFL)
64 DSHR=2D0*SQRT(DSH)
65 DSHZ=DSH/DBLE(ZS(JT))
66 XE=MAX(XE0,XB*(1./(1.-PARP(66))-1.))
67 IF(XB+XE.GE.0.999) THEN
68 Q2B=0.
69 GOTO 220
70 ENDIF
71
72C...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)
78 ELSE
79 Q2B=Q2S(JT)
80 TEVB=TEVS(JT)
81 ENDIF
82 ALSDUM=ULALPS_HIJING(PARP(63)*Q2B)
83 TEVB=TEVB+2.*LOG(ALAM(JT)/PARU(117))
84 TEVBSV=TEVB
85 ALAM(JT)=PARU(117)
86 B0=(33.-2.*MSTU(118))/6.
87
88C...Calculate Altarelli-Parisi and structure function weights.
89 DO 140 KFL=-6,6
90 WTAP(KFL)=0.
91 140 WTSF(KFL)=0.
92 IF(KFLB.EQ.21) THEN
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
97 ELSE
98 WTAP(0)=0.5*XB*(1./(XB+XE)-1.)
99 WTAP(KFLB)=8.*LOG((1.-XB)*(XB+XE)/XE)/3.
100 ENDIF
101 160 WTSUM=0.
102 IF(KFLB.NE.21) XFBO=XFB(KFLB)
103 IF(KFLB.EQ.21) XFBO=XFB(0)
104C***************************************************************
105C**********ERROR HAS OCCURED HERE
106 IF(XFBO.EQ.0.0) THEN
107 WRITE(MSTU(11),1000)
108 WRITE(MSTU(11),1001) KFLB,XFB(KFLB)
109 XFBO=0.00001
110 ENDIF
111C****************************************************************
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)
116
117C...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))
122 ELSE
123 TEVB=TEVB*EXP(MAX(-100.,LOG(RLU_HIJING(0))*B0/(5.*WTSUM)))
124 ENDIF
125 190 Q2REF=ALAM(JT)**2*EXP(TEVB)
126 Q2B=Q2REF/PARP(63)
127
128C...Evolution ended or select flavour for branching parton.
129 IF(Q2B.LT.PARP(62)**2) THEN
130 Q2B=0.
131 ELSE
132 WTRAN=RLU_HIJING(0)*WTSUM
133 KFLA=-MSTP(54)-1
134 200 KFLA=KFLA+1
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
138
139C...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))
142 WTZ=(1.-Z*(1.-Z))**2
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.))
148 WTZ=1.-2.*Z*(1.-Z)
149 ELSE
150 Z=1.-(1.-XB)*(XE/((XB+XE)*(1.-XB)))**RLU_HIJING(0)
151 WTZ=0.5*(1.+Z**2)
152 ENDIF
153
154C...Option with resummation of soft gluon emission as effective z shift.
155 IF(MSTP(65).GE.1) THEN
156 RSOFT=6.
157 IF(KFLB.NE.21) RSOFT=8./3.
158 Z=Z*(TEVB/TEVS(JT))**(RSOFT*XE/((XB+XE)*B0))
159 IF(Z.LE.XB) GOTO 180
160 ENDIF
161
162C...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.
168 ENDIF
169
170C...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
174 ENDIF
175
176C...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
182 TEVB=TEVBSV
183 WTAP(KFLB)=0.
184 GOTO 160
185 ELSEIF(TEVBSV-TEVB.GT.0.2) THEN
186 TEVB=0.5*(TEVBSV+TEVB)
187 GOTO 190
188 ELSE
189 XFBN=1E-10
190 ENDIF
191 ENDIF
192 DO 210 KFL=-MSTP(54),MSTP(54)
193 210 XFB(KFL)=XFN(KFL)
194 XA=XB/Z
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
202 ENDIF
203
204C...Define two hard scatterers in their CM-frame.
205 220 IF(N.EQ.NS+2) THEN
206 DQ2(JT)=Q2B
207 DPLCM=SQRT((DSH+DQ2(1)+DQ2(2))**2-4D0*DQ2(1)*DQ2(2))/DSHR
208 DO 240 JR=1,2
209 I=NS+JR
210 IF(JR.EQ.1) IPO=IPUS1
211 IF(JR.EQ.2) IPO=IPUS2
212 DO 230 J=1,5
213 K(I,J)=0
214 P(I,J)=0.
215 230 V(I,J)=0.
216 K(I,1)=14
217 K(I,2)=KFLS(JR+2)
218 K(I,4)=IPO
219 K(I,5)=IPO
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)))
223 K(IPO,1)=14
224 K(IPO,3)=I
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
227
228C...Find maximum allowed mass of timelike parton.
229 ELSEIF(N.GT.NS+2) THEN
230 JR=3-JT
231 DQ2(3)=Q2B
232 DPC(1)=P(IS(1),4)
233 DPC(2)=P(IS(2),4)
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))
239 IKIN=0
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)
246
247C...Generate timelike parton shower (if required).
248 IT=N
249 DO 250 J=1,5
250 K(IT,J)=0
251 P(IT,J)=0.
252 250 V(IT,J)=0.
253 K(IT,1)=3
254 K(IT,2)=21
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
263 Q2TIM=DMSMA
264 ELSEIF(MSTP(63).EQ.2) THEN
265 Q2TIM=MIN(SNGL(DMSMA),PARP(71)*Q2S(JT))
266 ELSE
267C'''Here remains to introduce angular ordering in first branching.
268 Q2TIM=DMSMA
269 ENDIF
270 CALL LUSHOW_HIJING(IT,0,SQRT(Q2TIM))
271 IF(N.GE.IT+1) P(IT,5)=P(IT+1,5)
272 ENDIF
273
274C...Reconstruct kinematics of branching: timelike parton shower.
275 DMS=P(IT,5)**2
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
285 IF(N.GE.IT+1) THEN
286 DPB(1)=SQRT(DPB(1)**2+DPT2)
287 DPB(2)=SQRT(DPB(1)**2+DMS)
288 DPB(3)=P(IT+1,3)
289 DPB(4)=SQRT(DPB(3)**2+DMS)
290 DBEZ=(DPB(4)*DPB(1)-DPB(3)*DPB(2))/(DPB(4)*DPB(2)-DPB(3)*
291 & DPB(1))
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)
295 ENDIF
296
297C...Reconstruct kinematics of branching: spacelike parton.
298 DO 260 J=1,5
299 K(N+1,J)=0
300 P(N+1,J)=0.
301 260 V(N+1,J)=0.
302 K(N+1,1)=14
303 K(N+1,2)=KFLB
304 P(N+1,1)=P(IT,1)
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)))
308
309C...Define colour flow of branching.
310 K(IS(JT),3)=N+1
311 K(IT,3)=N+1
312 ID1=IT
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)
317 ID2=IT+IS(JT)-ID1
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)
324 N=N+1
325
326C...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)
329 $ +P(IS(JR),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)
333 ENDIF
334
335C...Save quantities, loop back.
336 IS(JT)=N
337 Q2S(JT)=Q2B
338 DQ2(JT)=Q2B
339 IF(MSTP(62).GE.3) THE2(JT)=THE2T
340 DSH=DSHZ
341 IF(Q2B.GE.(0.5*PARP(62))**2) THEN
342 KFLS(JT+2)=KFLS(JT)
343 KFLS(JT)=KFLA
344 XS(JT)=XA
345 ZS(JT)=Z
346 DO 270 KFL=-6,6
347 270 XFS(JT,KFL)=XFA(KFL)
348 TEVS(JT)=TEVB
349 ELSE
350 IF(JT.EQ.1) IPU1=N
351 IF(JT.EQ.2) IPU2=N
352 ENDIF
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
358 ENDIF
359 IF(MAX(Q2S(1),Q2S(2)).GE.(0.5*PARP(62))**2.OR.N.LE.NS+1) GOTO 120
360
361C...Boost hard scattering partons to frame of shower initiators.
362 DO 280 J=1,3
363 280 ROBO(J+2)=(P(NS+1,J)+P(NS+2,J))/(P(NS+1,4)+P(NS+2,4))
364 DO 290 J=1,5
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
372 ENDIF
373 CALL LUDBRB_HIJING(N+2,N+2,0.,0.,-DBLE(ROBO(3)),-DBLE(ROBO(4)),
374 &-DBLE(ROBO(5)))
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)))
379
380C...Store user information. Reset Lambda value.
381 K(IPU1,3)=MINT(83)+3
382 K(IPU2,3)=MINT(83)+4
383 DO 300 JT=1,2
384 MINT(12+JT)=KFLS(JT)
385 300 VINT(140+JT)=XS(JT)
386 PARU(111)=ALAMS
387 1000 FORMAT(5X,'structure function has a zero point here')
388 1001 FORMAT(5X,'xf(x,i=',I5,')=',F10.5)
389
390 RETURN
391 END