New functions and constructors added and some other fixes.
[u/mrichter/AliRoot.git] / PYTHIA / jetset / lukfdi.F
1  
2 C********************************************************************* 
3  
4       SUBROUTINE LUKFDI(KFL1,KFL2,KFL3,KF) 
5  
6 C...Purpose: to generate a new flavour pair and combine off a hadron. 
7       COMMON/LUDAT1/MSTU(200),PARU(200),MSTJ(200),PARJ(200) 
8       COMMON/LUDAT2/KCHG(500,3),PMAS(500,4),PARF(2000),VCKM(4,4) 
9       SAVE /LUDAT1/,/LUDAT2/ 
10  
11 C...Default flavour values. Input consistency checks. 
12       KF1A=IABS(KFL1) 
13       KF2A=IABS(KFL2) 
14       KFL3=0 
15       KF=0 
16       IF(KF1A.EQ.0) RETURN 
17       IF(KF2A.NE.0) THEN 
18         IF(KF1A.LE.10.AND.KF2A.LE.10.AND.KFL1*KFL2.GT.0) RETURN 
19         IF(KF1A.GT.10.AND.KF2A.GT.10) RETURN 
20         IF((KF1A.GT.10.OR.KF2A.GT.10).AND.KFL1*KFL2.LT.0) RETURN 
21       ENDIF 
22  
23 C...Check if tabulated flavour probabilities are to be used. 
24       IF(MSTJ(15).EQ.1) THEN 
25         KTAB1=-1 
26         IF(KF1A.GE.1.AND.KF1A.LE.6) KTAB1=KF1A 
27         KFL1A=MOD(KF1A/1000,10) 
28         KFL1B=MOD(KF1A/100,10) 
29         KFL1S=MOD(KF1A,10) 
30         IF(KFL1A.GE.1.AND.KFL1A.LE.4.AND.KFL1B.GE.1.AND.KFL1B.LE.4) 
31      &  KTAB1=6+KFL1A*(KFL1A-2)+2*KFL1B+(KFL1S-1)/2 
32         IF(KFL1A.GE.1.AND.KFL1A.LE.4.AND.KFL1A.EQ.KFL1B) KTAB1=KTAB1-1 
33         IF(KF1A.GE.1.AND.KF1A.LE.6) KFL1A=KF1A 
34         KTAB2=0 
35         IF(KF2A.NE.0) THEN 
36           KTAB2=-1 
37           IF(KF2A.GE.1.AND.KF2A.LE.6) KTAB2=KF2A 
38           KFL2A=MOD(KF2A/1000,10) 
39           KFL2B=MOD(KF2A/100,10) 
40           KFL2S=MOD(KF2A,10) 
41           IF(KFL2A.GE.1.AND.KFL2A.LE.4.AND.KFL2B.GE.1.AND.KFL2B.LE.4) 
42      &    KTAB2=6+KFL2A*(KFL2A-2)+2*KFL2B+(KFL2S-1)/2 
43           IF(KFL2A.GE.1.AND.KFL2A.LE.4.AND.KFL2A.EQ.KFL2B) KTAB2=KTAB2-1 
44         ENDIF 
45         IF(KTAB1.GE.0.AND.KTAB2.GE.0) GOTO 150 
46       ENDIF 
47  
48 C...Parameters and breaking diquark parameter combinations. 
49   100 PAR2=PARJ(2) 
50       PAR3=PARJ(3) 
51       PAR4=3.*PARJ(4) 
52       IF(MSTJ(12).GE.2) THEN 
53         PAR3M=SQRT(PARJ(3)) 
54         PAR4M=1./(3.*SQRT(PARJ(4))) 
55         PARDM=PARJ(7)/(PARJ(7)+PAR3M*PARJ(6)) 
56         PARS0=PARJ(5)*(2.+(1.+PAR2*PAR3M*PARJ(7))*(1.+PAR4M)) 
57         PARS1=PARJ(7)*PARS0/(2.*PAR3M)+PARJ(5)*(PARJ(6)*(1.+PAR4M)+ 
58      &  PAR2*PAR3M*PARJ(6)*PARJ(7)) 
59         PARS2=PARJ(5)*2.*PARJ(6)*PARJ(7)*(PAR2*PARJ(7)+(1.+PAR4M)/PAR3M) 
60         PARSM=MAX(PARS0,PARS1,PARS2) 
61         PAR4=PAR4*(1.+PARSM)/(1.+PARSM/(3.*PAR4M)) 
62       ENDIF 
63  
64 C...Choice of whether to generate meson or baryon. 
65   110 MBARY=0 
66       KFDA=0 
67       IF(KF1A.LE.10) THEN 
68         IF(KF2A.EQ.0.AND.MSTJ(12).GE.1.AND.(1.+PARJ(1))*RLU(0).GT.1.) 
69      &  MBARY=1 
70         IF(KF2A.GT.10) MBARY=2 
71         IF(KF2A.GT.10.AND.KF2A.LE.10000) KFDA=KF2A 
72       ELSE 
73         MBARY=2 
74         IF(KF1A.LE.10000) KFDA=KF1A 
75       ENDIF 
76  
77 C...Possibility of process diquark -> meson + new diquark. 
78       IF(KFDA.NE.0.AND.MSTJ(12).GE.2) THEN 
79         KFLDA=MOD(KFDA/1000,10) 
80         KFLDB=MOD(KFDA/100,10) 
81         KFLDS=MOD(KFDA,10) 
82         WTDQ=PARS0 
83         IF(MAX(KFLDA,KFLDB).EQ.3) WTDQ=PARS1 
84         IF(MIN(KFLDA,KFLDB).EQ.3) WTDQ=PARS2 
85         IF(KFLDS.EQ.1) WTDQ=WTDQ/(3.*PAR4M) 
86         IF((1.+WTDQ)*RLU(0).GT.1.) MBARY=-1 
87         IF(MBARY.EQ.-1.AND.KF2A.NE.0) RETURN 
88       ENDIF 
89  
90 C...Flavour for meson, possibly with new flavour. 
91       IF(MBARY.LE.0) THEN 
92         KFS=ISIGN(1,KFL1) 
93         IF(MBARY.EQ.0) THEN 
94           IF(KF2A.EQ.0) KFL3=ISIGN(1+INT((2.+PAR2)*RLU(0)),-KFL1) 
95           KFLA=MAX(KF1A,KF2A+IABS(KFL3)) 
96           KFLB=MIN(KF1A,KF2A+IABS(KFL3)) 
97           IF(KFLA.NE.KF1A) KFS=-KFS 
98  
99 C...Splitting of diquark into meson plus new diquark. 
100         ELSE 
101           KFL1A=MOD(KF1A/1000,10) 
102           KFL1B=MOD(KF1A/100,10) 
103   120     KFL1D=KFL1A+INT(RLU(0)+0.5)*(KFL1B-KFL1A) 
104           KFL1E=KFL1A+KFL1B-KFL1D 
105           IF((KFL1D.EQ.3.AND.RLU(0).GT.PARDM).OR.(KFL1E.EQ.3.AND. 
106      &    RLU(0).LT.PARDM)) THEN 
107             KFL1D=KFL1A+KFL1B-KFL1D 
108             KFL1E=KFL1A+KFL1B-KFL1E 
109           ENDIF 
110           KFL3A=1+INT((2.+PAR2*PAR3M*PARJ(7))*RLU(0)) 
111           IF((KFL1E.NE.KFL3A.AND.RLU(0).GT.(1.+PAR4M)/MAX(2.,1.+PAR4M)) 
112      &    .OR.(KFL1E.EQ.KFL3A.AND.RLU(0).GT.2./MAX(2.,1.+PAR4M))) 
113      &    GOTO 120 
114           KFLDS=3 
115           IF(KFL1E.NE.KFL3A) KFLDS=2*INT(RLU(0)+1./(1.+PAR4M))+1 
116           KFL3=ISIGN(10000+1000*MAX(KFL1E,KFL3A)+100*MIN(KFL1E,KFL3A)+ 
117      &    KFLDS,-KFL1) 
118           KFLA=MAX(KFL1D,KFL3A) 
119           KFLB=MIN(KFL1D,KFL3A) 
120           IF(KFLA.NE.KFL1D) KFS=-KFS 
121         ENDIF 
122  
123 C...Form meson, with spin and flavour mixing for diagonal states. 
124         IF(KFLA.LE.2) KMUL=INT(PARJ(11)+RLU(0)) 
125         IF(KFLA.EQ.3) KMUL=INT(PARJ(12)+RLU(0)) 
126         IF(KFLA.GE.4) KMUL=INT(PARJ(13)+RLU(0)) 
127         IF(KMUL.EQ.0.AND.PARJ(14).GT.0.) THEN 
128           IF(RLU(0).LT.PARJ(14)) KMUL=2 
129         ELSEIF(KMUL.EQ.1.AND.PARJ(15)+PARJ(16)+PARJ(17).GT.0.) THEN 
130           RMUL=RLU(0) 
131           IF(RMUL.LT.PARJ(15)) KMUL=3 
132           IF(KMUL.EQ.1.AND.RMUL.LT.PARJ(15)+PARJ(16)) KMUL=4 
133           IF(KMUL.EQ.1.AND.RMUL.LT.PARJ(15)+PARJ(16)+PARJ(17)) KMUL=5 
134         ENDIF 
135         KFLS=3 
136         IF(KMUL.EQ.0.OR.KMUL.EQ.3) KFLS=1 
137         IF(KMUL.EQ.5) KFLS=5 
138         IF(KFLA.NE.KFLB) THEN 
139           KF=(100*KFLA+10*KFLB+KFLS)*KFS*(-1)**KFLA 
140         ELSE 
141           RMIX=RLU(0) 
142           IMIX=2*KFLA+10*KMUL 
143           IF(KFLA.LE.3) KF=110*(1+INT(RMIX+PARF(IMIX-1))+ 
144      &    INT(RMIX+PARF(IMIX)))+KFLS 
145           IF(KFLA.GE.4) KF=110*KFLA+KFLS 
146         ENDIF 
147         IF(KMUL.EQ.2.OR.KMUL.EQ.3) KF=KF+ISIGN(10000,KF) 
148         IF(KMUL.EQ.4) KF=KF+ISIGN(20000,KF) 
149  
150 C...Optional extra suppression of eta and eta'. 
151         IF(KF.EQ.221) THEN 
152           IF(RLU(0).GT.PARJ(25)) GOTO 110 
153         ELSEIF(KF.EQ.331) THEN 
154           IF(RLU(0).GT.PARJ(26)) GOTO 110 
155         ENDIF 
156  
157 C...Generate diquark flavour. 
158       ELSE 
159   130   IF(KF1A.LE.10.AND.KF2A.EQ.0) THEN 
160           KFLA=KF1A 
161   140     KFLB=1+INT((2.+PAR2*PAR3)*RLU(0)) 
162           KFLC=1+INT((2.+PAR2*PAR3)*RLU(0)) 
163           KFLDS=1 
164           IF(KFLB.GE.KFLC) KFLDS=3 
165           IF(KFLDS.EQ.1.AND.PAR4*RLU(0).GT.1.) GOTO 140 
166           IF(KFLDS.EQ.3.AND.PAR4.LT.RLU(0)) GOTO 140 
167           KFL3=ISIGN(1000*MAX(KFLB,KFLC)+100*MIN(KFLB,KFLC)+KFLDS,KFL1) 
168  
169 C...Take diquark flavour from input. 
170         ELSEIF(KF1A.LE.10) THEN 
171           KFLA=KF1A 
172           KFLB=MOD(KF2A/1000,10) 
173           KFLC=MOD(KF2A/100,10) 
174           KFLDS=MOD(KF2A,10) 
175  
176 C...Generate (or take from input) quark to go with diquark. 
177         ELSE 
178           IF(KF2A.EQ.0) KFL3=ISIGN(1+INT((2.+PAR2)*RLU(0)),KFL1) 
179           KFLA=KF2A+IABS(KFL3) 
180           KFLB=MOD(KF1A/1000,10) 
181           KFLC=MOD(KF1A/100,10) 
182           KFLDS=MOD(KF1A,10) 
183         ENDIF 
184  
185 C...SU(6) factors for formation of baryon. Try again if fails. 
186         KBARY=KFLDS 
187         IF(KFLDS.EQ.3.AND.KFLB.NE.KFLC) KBARY=5 
188         IF(KFLA.NE.KFLB.AND.KFLA.NE.KFLC) KBARY=KBARY+1 
189         WT=PARF(60+KBARY)+PARJ(18)*PARF(70+KBARY) 
190         IF(MBARY.EQ.1.AND.MSTJ(12).GE.2) THEN 
191           WTDQ=PARS0 
192           IF(MAX(KFLB,KFLC).EQ.3) WTDQ=PARS1 
193           IF(MIN(KFLB,KFLC).EQ.3) WTDQ=PARS2 
194           IF(KFLDS.EQ.1) WTDQ=WTDQ/(3.*PAR4M) 
195           IF(KFLDS.EQ.1) WT=WT*(1.+WTDQ)/(1.+PARSM/(3.*PAR4M)) 
196           IF(KFLDS.EQ.3) WT=WT*(1.+WTDQ)/(1.+PARSM) 
197         ENDIF 
198         IF(KF2A.EQ.0.AND.WT.LT.RLU(0)) GOTO 130 
199  
200 C...Form baryon. Distinguish Lambda- and Sigmalike baryons. 
201         KFLD=MAX(KFLA,KFLB,KFLC) 
202         KFLF=MIN(KFLA,KFLB,KFLC) 
203         KFLE=KFLA+KFLB+KFLC-KFLD-KFLF 
204         KFLS=2 
205         IF((PARF(60+KBARY)+PARJ(18)*PARF(70+KBARY))*RLU(0).GT. 
206      &  PARF(60+KBARY)) KFLS=4 
207         KFLL=0 
208         IF(KFLS.EQ.2.AND.KFLD.GT.KFLE.AND.KFLE.GT.KFLF) THEN 
209           IF(KFLDS.EQ.1.AND.KFLA.EQ.KFLD) KFLL=1 
210           IF(KFLDS.EQ.1.AND.KFLA.NE.KFLD) KFLL=INT(0.25+RLU(0)) 
211           IF(KFLDS.EQ.3.AND.KFLA.NE.KFLD) KFLL=INT(0.75+RLU(0)) 
212         ENDIF 
213         IF(KFLL.EQ.0) KF=ISIGN(1000*KFLD+100*KFLE+10*KFLF+KFLS,KFL1) 
214         IF(KFLL.EQ.1) KF=ISIGN(1000*KFLD+100*KFLF+10*KFLE+KFLS,KFL1) 
215       ENDIF 
216       RETURN 
217  
218 C...Use tabulated probabilities to select new flavour and hadron. 
219   150 IF(KTAB2.EQ.0.AND.MSTJ(12).LE.0) THEN 
220         KT3L=1 
221         KT3U=6 
222       ELSEIF(KTAB2.EQ.0.AND.KTAB1.GE.7.AND.MSTJ(12).LE.1) THEN 
223         KT3L=1 
224         KT3U=6 
225       ELSEIF(KTAB2.EQ.0) THEN 
226         KT3L=1 
227         KT3U=22 
228       ELSE 
229         KT3L=KTAB2 
230         KT3U=KTAB2 
231       ENDIF 
232       RFL=0. 
233       DO 170 KTS=0,2 
234       DO 160 KT3=KT3L,KT3U 
235       RFL=RFL+PARF(120+80*KTAB1+25*KTS+KT3) 
236   160 CONTINUE 
237   170 CONTINUE 
238       RFL=RLU(0)*RFL 
239       DO 190 KTS=0,2 
240       KTABS=KTS 
241       DO 180 KT3=KT3L,KT3U 
242       KTAB3=KT3 
243       RFL=RFL-PARF(120+80*KTAB1+25*KTS+KT3) 
244       IF(RFL.LE.0.) GOTO 200 
245   180 CONTINUE 
246   190 CONTINUE 
247   200 CONTINUE 
248  
249 C...Reconstruct flavour of produced quark/diquark. 
250       IF(KTAB3.LE.6) THEN 
251         KFL3A=KTAB3 
252         KFL3B=0 
253         KFL3=ISIGN(KFL3A,KFL1*(2*KTAB1-13)) 
254       ELSE 
255         KFL3A=1 
256         IF(KTAB3.GE.8) KFL3A=2 
257         IF(KTAB3.GE.11) KFL3A=3 
258         IF(KTAB3.GE.16) KFL3A=4 
259         KFL3B=(KTAB3-6-KFL3A*(KFL3A-2))/2 
260         KFL3=1000*KFL3A+100*KFL3B+1 
261         IF(KFL3A.EQ.KFL3B.OR.KTAB3.NE.6+KFL3A*(KFL3A-2)+2*KFL3B) KFL3= 
262      &  KFL3+2 
263         KFL3=ISIGN(KFL3,KFL1*(13-2*KTAB1)) 
264       ENDIF 
265  
266 C...Reconstruct meson code. 
267       IF(KFL3A.EQ.KFL1A.AND.KFL3B.EQ.KFL1B.AND.(KFL3A.LE.3.OR. 
268      &KFL3B.NE.0)) THEN 
269         RFL=RLU(0)*(PARF(143+80*KTAB1+25*KTABS)+PARF(144+80*KTAB1+ 
270      &  25*KTABS)+PARF(145+80*KTAB1+25*KTABS)) 
271         KF=110+2*KTABS+1 
272         IF(RFL.GT.PARF(143+80*KTAB1+25*KTABS)) KF=220+2*KTABS+1 
273         IF(RFL.GT.PARF(143+80*KTAB1+25*KTABS)+PARF(144+80*KTAB1+ 
274      &  25*KTABS)) KF=330+2*KTABS+1 
275       ELSEIF(KTAB1.LE.6.AND.KTAB3.LE.6) THEN 
276         KFLA=MAX(KTAB1,KTAB3) 
277         KFLB=MIN(KTAB1,KTAB3) 
278         KFS=ISIGN(1,KFL1) 
279         IF(KFLA.NE.KF1A) KFS=-KFS 
280         KF=(100*KFLA+10*KFLB+2*KTABS+1)*KFS*(-1)**KFLA 
281       ELSEIF(KTAB1.GE.7.AND.KTAB3.GE.7) THEN 
282         KFS=ISIGN(1,KFL1) 
283         IF(KFL1A.EQ.KFL3A) THEN 
284           KFLA=MAX(KFL1B,KFL3B) 
285           KFLB=MIN(KFL1B,KFL3B) 
286           IF(KFLA.NE.KFL1B) KFS=-KFS 
287         ELSEIF(KFL1A.EQ.KFL3B) THEN 
288           KFLA=KFL3A 
289           KFLB=KFL1B 
290           KFS=-KFS 
291         ELSEIF(KFL1B.EQ.KFL3A) THEN 
292           KFLA=KFL1A 
293           KFLB=KFL3B 
294         ELSEIF(KFL1B.EQ.KFL3B) THEN 
295           KFLA=MAX(KFL1A,KFL3A) 
296           KFLB=MIN(KFL1A,KFL3A) 
297           IF(KFLA.NE.KFL1A) KFS=-KFS 
298         ELSE 
299           CALL LUERRM(2,'(LUKFDI:) no matching flavours for qq -> qq') 
300           GOTO 100 
301         ENDIF 
302         KF=(100*KFLA+10*KFLB+2*KTABS+1)*KFS*(-1)**KFLA 
303  
304 C...Reconstruct baryon code. 
305       ELSE 
306         IF(KTAB1.GE.7) THEN 
307           KFLA=KFL3A 
308           KFLB=KFL1A 
309           KFLC=KFL1B 
310         ELSE 
311           KFLA=KFL1A 
312           KFLB=KFL3A 
313           KFLC=KFL3B 
314         ENDIF 
315         KFLD=MAX(KFLA,KFLB,KFLC) 
316         KFLF=MIN(KFLA,KFLB,KFLC) 
317         KFLE=KFLA+KFLB+KFLC-KFLD-KFLF 
318         IF(KTABS.EQ.0) KF=ISIGN(1000*KFLD+100*KFLF+10*KFLE+2,KFL1) 
319         IF(KTABS.GE.1) KF=ISIGN(1000*KFLD+100*KFLE+10*KFLF+2*KTABS,KFL1) 
320       ENDIF 
321  
322 C...Check that constructed flavour code is an allowed one. 
323       IF(KFL2.NE.0) KFL3=0 
324       KC=LUCOMP(KF) 
325       IF(KC.EQ.0) THEN 
326         CALL LUERRM(2,'(LUKFDI:) user-defined flavour probabilities '// 
327      &  'failed') 
328         GOTO 100 
329       ENDIF 
330  
331       RETURN 
332       END