]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PYTHIA/jetset/lukfdi.F
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Merging the VirtualMC branch to the main development branch (HEAD)
[u/mrichter/AliRoot.git] / PYTHIA / jetset / lukfdi.F
CommitLineData
fe4da5cc 1
2C*********************************************************************
3
4 SUBROUTINE LUKFDI(KFL1,KFL2,KFL3,KF)
5
6C...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
11C...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
23C...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
48C...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
64C...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
77C...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
90C...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
99C...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
123C...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
150C...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
157C...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
169C...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
176C...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
185C...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
200C...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
218C...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
249C...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
266C...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
304C...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
322C...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