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