]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PHOS/shaker/luprep.f
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Corrections thanks to A.Angelis
[u/mrichter/AliRoot.git] / PHOS / shaker / luprep.f
CommitLineData
fe4da5cc 1*CMZ : 17/07/98 15.44.31 by Federico Carminati
2*-- Author :
3C*********************************************************************
4
5 SUBROUTINE LUPREP(IP)
6
7C...Purpose: to rearrange partons along strings, to allow small systems
8C...to collapse into one or two particles and to check flavours.
9 IMPLICIT DOUBLE PRECISION(D)
10*KEEP,LUJETS.
11 COMMON /LUJETS/ N,K(200000,5),P(200000,5),V(200000,5)
12 SAVE /LUJETS/
13*KEEP,LUDAT1.
14 COMMON /LUDAT1/ MSTU(200),PARU(200),MSTJ(200),PARJ(200)
15 SAVE /LUDAT1/
16*KEEP,LUDAT2.
17 COMMON /LUDAT2/ KCHG(500,3),PMAS(500,4),PARF(2000),VCKM(4,4)
18 SAVE /LUDAT2/
19*KEEP,LUDAT3.
20 COMMON /LUDAT3/ MDCY(500,3),MDME(2000,2),BRAT(2000),KFDP(2000,5)
21 SAVE /LUDAT3/
22*KEND.
23 DIMENSION DPS(5),DPC(5),UE(3)
24
25C...Rearrange parton shower product listing along strings: begin loop.
26 I1=N
27 DO 130 MQGST=1,2
28 DO 120 I=MAX(1,IP),N
29 IF(K(I,1).NE.3) GOTO 120
30 KC=LUCOMP(K(I,2))
31 IF(KC.EQ.0) GOTO 120
32 KQ=KCHG(KC,2)
33 IF(KQ.EQ.0.OR.(MQGST.EQ.1.AND.KQ.EQ.2)) GOTO 120
34
35C...Pick up loose string end.
36 KCS=4
37 IF(KQ*ISIGN(1,K(I,2)).LT.0) KCS=5
38 IA=I
39 NSTP=0
40 100 NSTP=NSTP+1
41 IF(NSTP.GT.4*N) THEN
42 CALL LUERRM(14,'(LUPREP:) caught in infinite loop')
43 RETURN
44 ENDIF
45
46C...Copy undecayed parton.
47 IF(K(IA,1).EQ.3) THEN
48 IF(I1.GE.MSTU(4)-MSTU(32)-5) THEN
49 CALL LUERRM(11,'(LUPREP:) no more memory left in LUJETS')
50 RETURN
51 ENDIF
52 I1=I1+1
53 K(I1,1)=2
54 IF(NSTP.GE.2.AND.IABS(K(IA,2)).NE.21) K(I1,1)=1
55 K(I1,2)=K(IA,2)
56 K(I1,3)=IA
57 K(I1,4)=0
58 K(I1,5)=0
59 DO 110 J=1,5
60 P(I1,J)=P(IA,J)
61 110 V(I1,J)=V(IA,J)
62 K(IA,1)=K(IA,1)+10
63 IF(K(I1,1).EQ.1) GOTO 120
64 ENDIF
65
66C...Go to next parton in colour space.
67 IB=IA
68 IF(MOD(K(IB,KCS)/MSTU(5)**2,2).EQ.0.AND.MOD(K(IB,KCS),MSTU(5)).
69 &NE.0) THEN
70 IA=MOD(K(IB,KCS),MSTU(5))
71 K(IB,KCS)=K(IB,KCS)+MSTU(5)**2
72 MREV=0
73 ELSE
74 IF(K(IB,KCS).GE.2*MSTU(5)**2.OR.MOD(K(IB,KCS)/MSTU(5),MSTU(5)).
75 & EQ.0) KCS=9-KCS
76 IA=MOD(K(IB,KCS)/MSTU(5),MSTU(5))
77 K(IB,KCS)=K(IB,KCS)+2*MSTU(5)**2
78 MREV=1
79 ENDIF
80 IF(IA.LE.0.OR.IA.GT.N) THEN
81 CALL LUERRM(12,'(LUPREP:) colour rearrangement failed')
82 RETURN
83 ENDIF
84 IF(MOD(K(IA,4)/MSTU(5),MSTU(5)).EQ.IB.OR.MOD(K(IA,5)/MSTU(5),
85 &MSTU(5)).EQ.IB) THEN
86 IF(MREV.EQ.1) KCS=9-KCS
87 IF(MOD(K(IA,KCS)/MSTU(5),MSTU(5)).NE.IB) KCS=9-KCS
88 K(IA,KCS)=K(IA,KCS)+2*MSTU(5)**2
89 ELSE
90 IF(MREV.EQ.0) KCS=9-KCS
91 IF(MOD(K(IA,KCS),MSTU(5)).NE.IB) KCS=9-KCS
92 K(IA,KCS)=K(IA,KCS)+MSTU(5)**2
93 ENDIF
94 IF(IA.NE.I) GOTO 100
95 K(I1,1)=1
96 120 CONTINUE
97 130 CONTINUE
98 N=I1
99
100C...Find lowest-mass colour singlet jet system, OK if above threshold.
101 IF(MSTJ(14).LE.0) GOTO 320
102 NS=N
103 140 NSIN=N-NS
104 PDM=1.+PARJ(32)
105 IC=0
106 DO 190 I=MAX(1,IP),NS
107 IF(K(I,1).NE.1.AND.K(I,1).NE.2) THEN
108 ELSEIF(K(I,1).EQ.2.AND.IC.EQ.0) THEN
109 NSIN=NSIN+1
110 IC=I
111 DO 150 J=1,4
112 150 DPS(J)=P(I,J)
113 MSTJ(93)=1
114 DPS(5)=ULMASS(K(I,2))
115 ELSEIF(K(I,1).EQ.2) THEN
116 DO 160 J=1,4
117 160 DPS(J)=DPS(J)+P(I,J)
118 ELSEIF(IC.NE.0.AND.KCHG(LUCOMP(K(I,2)),2).NE.0) THEN
119 DO 170 J=1,4
120 170 DPS(J)=DPS(J)+P(I,J)
121 MSTJ(93)=1
122 DPS(5)=DPS(5)+ULMASS(K(I,2))
123 PD=SQRT(MAX(0D0,DPS(4)**2-DPS(1)**2-DPS(2)**2-DPS(3)**2))-DPS(5)
124 IF(PD.LT.PDM) THEN
125 PDM=PD
126 DO 180 J=1,5
127 180 DPC(J)=DPS(J)
128 IC1=IC
129 IC2=I
130 ENDIF
131 IC=0
132 ELSE
133 NSIN=NSIN+1
134 ENDIF
135 190 CONTINUE
136 IF(PDM.GE.PARJ(32)) GOTO 320
137
138C...Fill small-mass system as cluster.
139 NSAV=N
140 PECM=SQRT(MAX(0D0,DPC(4)**2-DPC(1)**2-DPC(2)**2-DPC(3)**2))
141 K(N+1,1)=11
142 K(N+1,2)=91
143 K(N+1,3)=IC1
144 K(N+1,4)=N+2
145 K(N+1,5)=N+3
146 P(N+1,1)=DPC(1)
147 P(N+1,2)=DPC(2)
148 P(N+1,3)=DPC(3)
149 P(N+1,4)=DPC(4)
150 P(N+1,5)=PECM
151
152C...Form two particles from flavours of lowest-mass system, if feasible.
153 K(N+2,1)=1
154 K(N+3,1)=1
155 IF(MSTU(16).NE.2) THEN
156 K(N+2,3)=N+1
157 K(N+3,3)=N+1
158 ELSE
159 K(N+2,3)=IC1
160 K(N+3,3)=IC2
161 ENDIF
162 K(N+2,4)=0
163 K(N+3,4)=0
164 K(N+2,5)=0
165 K(N+3,5)=0
166 IF(IABS(K(IC1,2)).NE.21) THEN
167 KC1=LUCOMP(K(IC1,2))
168 KC2=LUCOMP(K(IC2,2))
169 IF(KC1.EQ.0.OR.KC2.EQ.0) GOTO 320
170 KQ1=KCHG(KC1,2)*ISIGN(1,K(IC1,2))
171 KQ2=KCHG(KC2,2)*ISIGN(1,K(IC2,2))
172 IF(KQ1+KQ2.NE.0) GOTO 320
173 200 CALL LUKFDI(K(IC1,2),0,KFLN,K(N+2,2))
174 CALL LUKFDI(K(IC2,2),-KFLN,KFLDMP,K(N+3,2))
175 IF(K(N+2,2).EQ.0.OR.K(N+3,2).EQ.0) GOTO 200
176 ELSE
177 IF(IABS(K(IC2,2)).NE.21) GOTO 320
178 210 CALL LUKFDI(1+INT((2.+PARJ(2))*RLU(0)),0,KFLN,KFDMP)
179 CALL LUKFDI(KFLN,0,KFLM,K(N+2,2))
180 CALL LUKFDI(-KFLN,-KFLM,KFLDMP,K(N+3,2))
181 IF(K(N+2,2).EQ.0.OR.K(N+3,2).EQ.0) GOTO 210
182 ENDIF
183 P(N+2,5)=ULMASS(K(N+2,2))
184 P(N+3,5)=ULMASS(K(N+3,2))
185 IF(P(N+2,5)+P(N+3,5)+PARJ(64).GE.PECM.AND.NSIN.EQ.1) GOTO 320
186 IF(P(N+2,5)+P(N+3,5)+PARJ(64).GE.PECM) GOTO 260
187
188C...Perform two-particle decay of jet system, if possible.
189 IF(PECM.GE.0.02*DPC(4)) THEN
190 PA=SQRT((PECM**2-(P(N+2,5)+P(N+3,5))**2)*(PECM**2-
191 & (P(N+2,5)-P(N+3,5))**2))/(2.*PECM)
192 UE(3)=2.*RLU(0)-1.
193 PHI=PARU(2)*RLU(0)
194 UE(1)=SQRT(1.-UE(3)**2)*COS(PHI)
195 UE(2)=SQRT(1.-UE(3)**2)*SIN(PHI)
196 DO 220 J=1,3
197 P(N+2,J)=PA*UE(J)
198 220 P(N+3,J)=-PA*UE(J)
199 P(N+2,4)=SQRT(PA**2+P(N+2,5)**2)
200 P(N+3,4)=SQRT(PA**2+P(N+3,5)**2)
201 MSTU(33)=1
202 CALL LUDBRB(N+2,N+3,0.,0.,DPC(1)/DPC(4),DPC(2)/DPC(4),
203 & DPC(3)/DPC(4))
204 ELSE
205 NP=0
206 DO 230 I=IC1,IC2
207 230 IF(K(I,1).EQ.1.OR.K(I,1).EQ.2) NP=NP+1
208 HA=P(IC1,4)*P(IC2,4)-P(IC1,1)*P(IC2,1)-P(IC1,2)*P(IC2,2)-
209 & P(IC1,3)*P(IC2,3)
210 IF(NP.GE.3.OR.HA.LE.1.25*P(IC1,5)*P(IC2,5)) GOTO 260
211 HD1=0.5*(P(N+2,5)**2-P(IC1,5)**2)
212 HD2=0.5*(P(N+3,5)**2-P(IC2,5)**2)
213 HR=SQRT(MAX(0.,((HA-HD1-HD2)**2-(P(N+2,5)*P(N+3,5))**2)/
214 & (HA**2-(P(IC1,5)*P(IC2,5))**2)))-1.
215 HC=P(IC1,5)**2+2.*HA+P(IC2,5)**2
216 HK1=((P(IC2,5)**2+HA)*HR+HD1-HD2)/HC
217 HK2=((P(IC1,5)**2+HA)*HR+HD2-HD1)/HC
218 DO 240 J=1,4
219 P(N+2,J)=(1.+HK1)*P(IC1,J)-HK2*P(IC2,J)
220 240 P(N+3,J)=(1.+HK2)*P(IC2,J)-HK1*P(IC1,J)
221 ENDIF
222 DO 250 J=1,4
223 V(N+1,J)=V(IC1,J)
224 V(N+2,J)=V(IC1,J)
225 250 V(N+3,J)=V(IC2,J)
226 V(N+1,5)=0.
227 V(N+2,5)=0.
228 V(N+3,5)=0.
229 N=N+3
230 GOTO 300
231
232C...Else form one particle from the flavours available, if possible.
233 260 K(N+1,5)=N+2
234 IF(IABS(K(IC1,2)).GT.100.AND.IABS(K(IC2,2)).GT.100) THEN
235 GOTO 320
236 ELSEIF(IABS(K(IC1,2)).NE.21) THEN
237 CALL LUKFDI(K(IC1,2),K(IC2,2),KFLDMP,K(N+2,2))
238 ELSE
239 KFLN=1+INT((2.+PARJ(2))*RLU(0))
240 CALL LUKFDI(KFLN,-KFLN,KFLDMP,K(N+2,2))
241 ENDIF
242 IF(K(N+2,2).EQ.0) GOTO 260
243 P(N+2,5)=ULMASS(K(N+2,2))
244
245C...Find parton/particle which combines to largest extra mass.
246 IR=0
247 HA=0.
248 HSM=0.
249 DO 280 MCOMB=1,3
250 IF(IR.NE.0) GOTO 280
251 DO 270 I=MAX(1,IP),N
252 IF(K(I,1).LE.0.OR.K(I,1).GT.10.OR.(I.GE.IC1.AND.I.LE.IC2.
253 &AND.K(I,1).GE.1.AND.K(I,1).LE.2)) GOTO 270
254 IF(MCOMB.EQ.1) KCI=LUCOMP(K(I,2))
255 IF(MCOMB.EQ.1.AND.KCI.EQ.0) GOTO 270
256 IF(MCOMB.EQ.1.AND.KCHG(KCI,2).EQ.0.AND.I.LE.NS) GOTO 270
257 IF(MCOMB.EQ.2.AND.IABS(K(I,2)).GT.10.AND.IABS(K(I,2)).LE.100)
258 &GOTO 270
259 HCR=DPC(4)*P(I,4)-DPC(1)*P(I,1)-DPC(2)*P(I,2)-DPC(3)*P(I,3)
260 HSR=2.*HCR+PECM**2-P(N+2,5)**2-2.*P(N+2,5)*P(I,5)
261 IF(HSR.GT.HSM) THEN
262 IR=I
263 HA=HCR
264 HSM=HSR
265 ENDIF
266 270 CONTINUE
267 280 CONTINUE
268
269C...Shuffle energy and momentum to put new particle on mass shell.
270 IF(IR.NE.0) THEN
271 HB=PECM**2+HA
272 HC=P(N+2,5)**2+HA
273 HD=P(IR,5)**2+HA
274 HK2=0.5*(HB*SQRT(MAX(0.,((HB+HC)**2-4.*(HB+HD)*P(N+2,5)**2)/
275 & (HA**2-(PECM*P(IR,5))**2)))-(HB+HC))/(HB+HD)
276 HK1=(0.5*(P(N+2,5)**2-PECM**2)+HD*HK2)/HB
277 DO 290 J=1,4
278 P(N+2,J)=(1.+HK1)*DPC(J)-HK2*P(IR,J)
279 P(IR,J)=(1.+HK2)*P(IR,J)-HK1*DPC(J)
280 V(N+1,J)=V(IC1,J)
281 290 V(N+2,J)=V(IC1,J)
282 V(N+1,5)=0.
283 V(N+2,5)=0.
284 N=N+2
285 ELSE
286 CALL LUERRM(3,'(LUPREP:) no match for collapsing cluster')
287 RETURN
288 ENDIF
289
290C...Mark collapsed system and store daughter pointers. Iterate.
291 300 DO 310 I=IC1,IC2
292 IF((K(I,1).EQ.1.OR.K(I,1).EQ.2).AND.KCHG(LUCOMP(K(I,2)),2).NE.0)
293 &THEN
294 K(I,1)=K(I,1)+10
295 IF(MSTU(16).NE.2) THEN
296 K(I,4)=NSAV+1
297 K(I,5)=NSAV+1
298 ELSE
299 K(I,4)=NSAV+2
300 K(I,5)=N
301 ENDIF
302 ENDIF
303 310 CONTINUE
304 IF(N.LT.MSTU(4)-MSTU(32)-5) GOTO 140
305
306C...Check flavours and invariant masses in parton systems.
307 320 NP=0
308 KFN=0
309 KQS=0
310 DO 330 J=1,5
311 330 DPS(J)=0.
312 DO 360 I=MAX(1,IP),N
313 IF(K(I,1).LE.0.OR.K(I,1).GT.10) GOTO 360
314 KC=LUCOMP(K(I,2))
315 IF(KC.EQ.0) GOTO 360
316 KQ=KCHG(KC,2)*ISIGN(1,K(I,2))
317 IF(KQ.EQ.0) GOTO 360
318 NP=NP+1
319 IF(KQ.NE.2) THEN
320 KFN=KFN+1
321 KQS=KQS+KQ
322 MSTJ(93)=1
323 DPS(5)=DPS(5)+ULMASS(K(I,2))
324 ENDIF
325 DO 340 J=1,4
326 340 DPS(J)=DPS(J)+P(I,J)
327 IF(K(I,1).EQ.1) THEN
328 IF(NP.NE.1.AND.(KFN.EQ.1.OR.KFN.GE.3.OR.KQS.NE.0)) CALL
329 & LUERRM(2,'(LUPREP:) unphysical flavour combination')
330 IF(NP.NE.1.AND.DPS(4)**2-DPS(1)**2-DPS(2)**2-DPS(3)**2.LT.
331 & (0.9*PARJ(32)+DPS(5))**2) CALL LUERRM(3,
332 & '(LUPREP:) too small mass in jet system')
333 NP=0
334 KFN=0
335 KQS=0
336 DO 350 J=1,5
337 350 DPS(J)=0.
338 ENDIF
339 360 CONTINUE
340
341 RETURN
342 END