Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
lukfdi.f
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file 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  SAVE /ludat1/
9  common/ludat2/kchg(500,3),pmas(500,4),parf(2000),vckm(4,4)
10  SAVE /ludat2/
11 
12 C...Default flavour values. Input consistency checks.
13  kf1a=iabs(kfl1)
14  kf2a=iabs(kfl2)
15  kfl3=0
16  kf=0
17  IF(kf1a.EQ.0) RETURN
18  IF(kf2a.NE.0) THEN
19  IF(kf1a.LE.10.AND.kf2a.LE.10.AND.kfl1*kfl2.GT.0) RETURN
20  IF(kf1a.GT.10.AND.kf2a.GT.10) RETURN
21  IF((kf1a.GT.10.OR.kf2a.GT.10).AND.kfl1*kfl2.LT.0) RETURN
22  ENDIF
23 
24 C...Check if tabulated flavour probabilities are to be used.
25  IF(mstj(15).EQ.1) THEN
26  ktab1=-1
27  IF(kf1a.GE.1.AND.kf1a.LE.6) ktab1=kf1a
28  kfl1a=mod(kf1a/1000,10)
29  kfl1b=mod(kf1a/100,10)
30  kfl1s=mod(kf1a,10)
31  IF(kfl1a.GE.1.AND.kfl1a.LE.4.AND.kfl1b.GE.1.AND.kfl1b.LE.4)
32  & ktab1=6+kfl1a*(kfl1a-2)+2*kfl1b+(kfl1s-1)/2
33  IF(kfl1a.GE.1.AND.kfl1a.LE.4.AND.kfl1a.EQ.kfl1b) ktab1=ktab1-1
34  IF(kf1a.GE.1.AND.kf1a.LE.6) kfl1a=kf1a
35  ktab2=0
36  IF(kf2a.NE.0) THEN
37  ktab2=-1
38  IF(kf2a.GE.1.AND.kf2a.LE.6) ktab2=kf2a
39  kfl2a=mod(kf2a/1000,10)
40  kfl2b=mod(kf2a/100,10)
41  kfl2s=mod(kf2a,10)
42  IF(kfl2a.GE.1.AND.kfl2a.LE.4.AND.kfl2b.GE.1.AND.kfl2b.LE.4)
43  & ktab2=6+kfl2a*(kfl2a-2)+2*kfl2b+(kfl2s-1)/2
44  IF(kfl2a.GE.1.AND.kfl2a.LE.4.AND.kfl2a.EQ.kfl2b) ktab2=ktab2-1
45  ENDIF
46  IF(ktab1.GE.0.AND.ktab2.GE.0) goto 140
47  ENDIF
48 
49 C...Parameters and breaking diquark parameter combinations.
50  100 par2=parj(2)
51  par3=parj(3)
52  par4=3.*parj(4)
53  IF(mstj(12).GE.2) THEN
54  par3m=sqrt(parj(3))
55  par4m=1./(3.*sqrt(parj(4)))
56  pardm=parj(7)/(parj(7)+par3m*parj(6))
57  pars0=parj(5)*(2.+(1.+par2*par3m*parj(7))*(1.+par4m))
58  pars1=parj(7)*pars0/(2.*par3m)+parj(5)*(parj(6)*(1.+par4m)+
59  & par2*par3m*parj(6)*parj(7))
60  pars2=parj(5)*2.*parj(6)*parj(7)*(par2*parj(7)+(1.+par4m)/par3m)
61  parsm=max(pars0,pars1,pars2)
62  par4=par4*(1.+parsm)/(1.+parsm/(3.*par4m))
63  ENDIF
64 
65 C...Choice of whether to generate meson or baryon.
66  mbary=0
67  kfda=0
68  IF(kf1a.LE.10) THEN
69  IF(kf2a.EQ.0.AND.mstj(12).GE.1.AND.(1.+parj(1))*rlu(0).GT.1.)
70  & mbary=1
71  IF(kf2a.GT.10) mbary=2
72  IF(kf2a.GT.10.AND.kf2a.LE.10000) kfda=kf2a
73  ELSE
74  mbary=2
75  IF(kf1a.LE.10000) kfda=kf1a
76  ENDIF
77 
78 C...Possibility of process diquark -> meson + new diquark.
79  IF(kfda.NE.0.AND.mstj(12).GE.2) THEN
80  kflda=mod(kfda/1000,10)
81  kfldb=mod(kfda/100,10)
82  kflds=mod(kfda,10)
83  wtdq=pars0
84  IF(max(kflda,kfldb).EQ.3) wtdq=pars1
85  IF(min(kflda,kfldb).EQ.3) wtdq=pars2
86  IF(kflds.EQ.1) wtdq=wtdq/(3.*par4m)
87  IF((1.+wtdq)*rlu(0).GT.1.) mbary=-1
88  IF(mbary.EQ.-1.AND.kf2a.NE.0) RETURN
89  ENDIF
90 
91 C...Flavour for meson, possibly with new flavour.
92  IF(mbary.LE.0) THEN
93  kfs=isign(1,kfl1)
94  IF(mbary.EQ.0) THEN
95  IF(kf2a.EQ.0) kfl3=isign(1+int((2.+par2)*rlu(0)),-kfl1)
96  kfla=max(kf1a,kf2a+iabs(kfl3))
97  kflb=min(kf1a,kf2a+iabs(kfl3))
98  IF(kfla.NE.kf1a) kfs=-kfs
99 
100 C...Splitting of diquark into meson plus new diquark.
101  ELSE
102  kfl1a=mod(kf1a/1000,10)
103  kfl1b=mod(kf1a/100,10)
104  110 kfl1d=kfl1a+int(rlu(0)+0.5)*(kfl1b-kfl1a)
105  kfl1e=kfl1a+kfl1b-kfl1d
106  IF((kfl1d.EQ.3.AND.rlu(0).GT.pardm).OR.(kfl1e.EQ.3.AND.
107  & rlu(0).LT.pardm)) THEN
108  kfl1d=kfl1a+kfl1b-kfl1d
109  kfl1e=kfl1a+kfl1b-kfl1e
110  ENDIF
111  kfl3a=1+int((2.+par2*par3m*parj(7))*rlu(0))
112  IF((kfl1e.NE.kfl3a.AND.rlu(0).GT.(1.+par4m)/max(2.,1.+par4m)).
113  & or.(kfl1e.EQ.kfl3a.AND.rlu(0).GT.2./max(2.,1.+par4m)))
114  & goto 110
115  kflds=3
116  IF(kfl1e.NE.kfl3a) kflds=2*int(rlu(0)+1./(1.+par4m))+1
117  kfl3=isign(10000+1000*max(kfl1e,kfl3a)+100*min(kfl1e,kfl3a)+
118  & kflds,-kfl1)
119  kfla=max(kfl1d,kfl3a)
120  kflb=min(kfl1d,kfl3a)
121  IF(kfla.NE.kfl1d) kfs=-kfs
122  ENDIF
123 
124 C...Form meson, with spin and flavour mixing for diagonal states.
125  IF(kfla.LE.2) kmul=int(parj(11)+rlu(0))
126  IF(kfla.EQ.3) kmul=int(parj(12)+rlu(0))
127  IF(kfla.GE.4) kmul=int(parj(13)+rlu(0))
128  IF(kmul.EQ.0.AND.parj(14).GT.0.) THEN
129  IF(rlu(0).LT.parj(14)) kmul=2
130  ELSEIF(kmul.EQ.1.AND.parj(15)+parj(16)+parj(17).GT.0.) THEN
131  rmul=rlu(0)
132  IF(rmul.LT.parj(15)) kmul=3
133  IF(kmul.EQ.1.AND.rmul.LT.parj(15)+parj(16)) kmul=4
134  IF(kmul.EQ.1.AND.rmul.LT.parj(15)+parj(16)+parj(17)) kmul=5
135  ENDIF
136  kfls=3
137  IF(kmul.EQ.0.OR.kmul.EQ.3) kfls=1
138  IF(kmul.EQ.5) kfls=5
139  IF(kfla.NE.kflb) THEN
140  kf=(100*kfla+10*kflb+kfls)*kfs*(-1)**kfla
141  ELSE
142  rmix=rlu(0)
143  imix=2*kfla+10*kmul
144  IF(kfla.LE.3) kf=110*(1+int(rmix+parf(imix-1))+
145  & int(rmix+parf(imix)))+kfls
146  IF(kfla.GE.4) kf=110*kfla+kfls
147  ENDIF
148  IF(kmul.EQ.2.OR.kmul.EQ.3) kf=kf+isign(10000,kf)
149  IF(kmul.EQ.4) kf=kf+isign(20000,kf)
150 
151 C...Generate diquark flavour.
152  ELSE
153  120 IF(kf1a.LE.10.AND.kf2a.EQ.0) THEN
154  kfla=kf1a
155  130 kflb=1+int((2.+par2*par3)*rlu(0))
156  kflc=1+int((2.+par2*par3)*rlu(0))
157  kflds=1
158  IF(kflb.GE.kflc) kflds=3
159  IF(kflds.EQ.1.AND.par4*rlu(0).GT.1.) goto 130
160  IF(kflds.EQ.3.AND.par4.LT.rlu(0)) goto 130
161  kfl3=isign(1000*max(kflb,kflc)+100*min(kflb,kflc)+kflds,kfl1)
162 
163 C...Take diquark flavour from input.
164  ELSEIF(kf1a.LE.10) THEN
165  kfla=kf1a
166  kflb=mod(kf2a/1000,10)
167  kflc=mod(kf2a/100,10)
168  kflds=mod(kf2a,10)
169 
170 C...Generate (or take from input) quark to go with diquark.
171  ELSE
172  IF(kf2a.EQ.0) kfl3=isign(1+int((2.+par2)*rlu(0)),kfl1)
173  kfla=kf2a+iabs(kfl3)
174  kflb=mod(kf1a/1000,10)
175  kflc=mod(kf1a/100,10)
176  kflds=mod(kf1a,10)
177  ENDIF
178 
179 C...SU(6) factors for formation of baryon. Try again if fails.
180  kbary=kflds
181  IF(kflds.EQ.3.AND.kflb.NE.kflc) kbary=5
182  IF(kfla.NE.kflb.AND.kfla.NE.kflc) kbary=kbary+1
183  wt=parf(60+kbary)+parj(18)*parf(70+kbary)
184  IF(mbary.EQ.1.AND.mstj(12).GE.2) THEN
185  wtdq=pars0
186  IF(max(kflb,kflc).EQ.3) wtdq=pars1
187  IF(min(kflb,kflc).EQ.3) wtdq=pars2
188  IF(kflds.EQ.1) wtdq=wtdq/(3.*par4m)
189  IF(kflds.EQ.1) wt=wt*(1.+wtdq)/(1.+parsm/(3.*par4m))
190  IF(kflds.EQ.3) wt=wt*(1.+wtdq)/(1.+parsm)
191  ENDIF
192  IF(kf2a.EQ.0.AND.wt.LT.rlu(0)) goto 120
193 
194 C...Form baryon. Distinguish Lambda- and Sigmalike baryons.
195  kfld=max(kfla,kflb,kflc)
196  kflf=min(kfla,kflb,kflc)
197  kfle=kfla+kflb+kflc-kfld-kflf
198  kfls=2
199  IF((parf(60+kbary)+parj(18)*parf(70+kbary))*rlu(0).GT.
200  & parf(60+kbary)) kfls=4
201  kfll=0
202  IF(kfls.EQ.2.AND.kfld.GT.kfle.AND.kfle.GT.kflf) THEN
203  IF(kflds.EQ.1.AND.kfla.EQ.kfld) kfll=1
204  IF(kflds.EQ.1.AND.kfla.NE.kfld) kfll=int(0.25+rlu(0))
205  IF(kflds.EQ.3.AND.kfla.NE.kfld) kfll=int(0.75+rlu(0))
206  ENDIF
207  IF(kfll.EQ.0) kf=isign(1000*kfld+100*kfle+10*kflf+kfls,kfl1)
208  IF(kfll.EQ.1) kf=isign(1000*kfld+100*kflf+10*kfle+kfls,kfl1)
209  ENDIF
210  RETURN
211 
212 C...Use tabulated probabilities to select new flavour and hadron.
213  140 IF(ktab2.EQ.0.AND.mstj(12).LE.0) THEN
214  kt3l=1
215  kt3u=6
216  ELSEIF(ktab2.EQ.0.AND.ktab1.GE.7.AND.mstj(12).LE.1) THEN
217  kt3l=1
218  kt3u=6
219  ELSEIF(ktab2.EQ.0) THEN
220  kt3l=1
221  kt3u=22
222  ELSE
223  kt3l=ktab2
224  kt3u=ktab2
225  ENDIF
226  rfl=0.
227  DO 150 kts=0,2
228  DO 150 kt3=kt3l,kt3u
229  rfl=rfl+parf(120+80*ktab1+25*kts+kt3)
230  150 CONTINUE
231  rfl=rlu(0)*rfl
232  DO 160 kts=0,2
233  ktabs=kts
234  DO 160 kt3=kt3l,kt3u
235  ktab3=kt3
236  rfl=rfl-parf(120+80*ktab1+25*kts+kt3)
237  160 IF(rfl.LE.0.) goto 170
238  170 CONTINUE
239 
240 C...Reconstruct flavour of produced quark/diquark.
241  IF(ktab3.LE.6) THEN
242  kfl3a=ktab3
243  kfl3b=0
244  kfl3=isign(kfl3a,kfl1*(2*ktab1-13))
245  ELSE
246  kfl3a=1
247  IF(ktab3.GE.8) kfl3a=2
248  IF(ktab3.GE.11) kfl3a=3
249  IF(ktab3.GE.16) kfl3a=4
250  kfl3b=(ktab3-6-kfl3a*(kfl3a-2))/2
251  kfl3=1000*kfl3a+100*kfl3b+1
252  IF(kfl3a.EQ.kfl3b.OR.ktab3.NE.6+kfl3a*(kfl3a-2)+2*kfl3b) kfl3=
253  & kfl3+2
254  kfl3=isign(kfl3,kfl1*(13-2*ktab1))
255  ENDIF
256 
257 C...Reconstruct meson code.
258  IF(kfl3a.EQ.kfl1a.AND.kfl3b.EQ.kfl1b.AND.(kfl3a.LE.3.OR.
259  &kfl3b.NE.0)) THEN
260  rfl=rlu(0)*(parf(143+80*ktab1+25*ktabs)+parf(144+80*ktab1+
261  & 25*ktabs)+parf(145+80*ktab1+25*ktabs))
262  kf=110+2*ktabs+1
263  IF(rfl.GT.parf(143+80*ktab1+25*ktabs)) kf=220+2*ktabs+1
264  IF(rfl.GT.parf(143+80*ktab1+25*ktabs)+parf(144+80*ktab1+
265  & 25*ktabs)) kf=330+2*ktabs+1
266  ELSEIF(ktab1.LE.6.AND.ktab3.LE.6) THEN
267  kfla=max(ktab1,ktab3)
268  kflb=min(ktab1,ktab3)
269  kfs=isign(1,kfl1)
270  IF(kfla.NE.kf1a) kfs=-kfs
271  kf=(100*kfla+10*kflb+2*ktabs+1)*kfs*(-1)**kfla
272  ELSEIF(ktab1.GE.7.AND.ktab3.GE.7) THEN
273  kfs=isign(1,kfl1)
274  IF(kfl1a.EQ.kfl3a) THEN
275  kfla=max(kfl1b,kfl3b)
276  kflb=min(kfl1b,kfl3b)
277  IF(kfla.NE.kfl1b) kfs=-kfs
278  ELSEIF(kfl1a.EQ.kfl3b) THEN
279  kfla=kfl3a
280  kflb=kfl1b
281  kfs=-kfs
282  ELSEIF(kfl1b.EQ.kfl3a) THEN
283  kfla=kfl1a
284  kflb=kfl3b
285  ELSEIF(kfl1b.EQ.kfl3b) THEN
286  kfla=max(kfl1a,kfl3a)
287  kflb=min(kfl1a,kfl3a)
288  IF(kfla.NE.kfl1a) kfs=-kfs
289  ELSE
290  CALL luerrm(2,'(LUKFDI:) no matching flavours for qq -> qq')
291  goto 100
292  ENDIF
293  kf=(100*kfla+10*kflb+2*ktabs+1)*kfs*(-1)**kfla
294 
295 C...Reconstruct baryon code.
296  ELSE
297  IF(ktab1.GE.7) THEN
298  kfla=kfl3a
299  kflb=kfl1a
300  kflc=kfl1b
301  ELSE
302  kfla=kfl1a
303  kflb=kfl3a
304  kflc=kfl3b
305  ENDIF
306  kfld=max(kfla,kflb,kflc)
307  kflf=min(kfla,kflb,kflc)
308  kfle=kfla+kflb+kflc-kfld-kflf
309  IF(ktabs.EQ.0) kf=isign(1000*kfld+100*kflf+10*kfle+2,kfl1)
310  IF(ktabs.GE.1) kf=isign(1000*kfld+100*kfle+10*kflf+2*ktabs,kfl1)
311  ENDIF
312 
313 C...Check that constructed flavour code is an allowed one.
314  IF(kfl2.NE.0) kfl3=0
315  kc=lucomp(kf)
316  IF(kc.EQ.0) THEN
317  CALL luerrm(2,'(LUKFDI:) user-defined flavour probabilities '//
318  & 'failed')
319  goto 100
320  ENDIF
321 
322  RETURN
323  END