]> git.uio.no Git - u/mrichter/AliRoot.git/blob - HIJING/hipyset1_35/lukfdi_hijing.F
Fixing the decoding of regional header bits.
[u/mrichter/AliRoot.git] / HIJING / hipyset1_35 / lukfdi_hijing.F
1 * $Id$
2     
3 C*********************************************************************  
4     
5       SUBROUTINE LUKFDI_HIJING(KFL1,KFL2,KFL3,KF)  
6     
7 C...Purpose: to generate a new flavour pair and combine off a hadron.   
8 #include "ludat1_hijing.inc"
9 #include "ludat2_hijing.inc"
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 140  
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       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_HIJING(0)
69      $        .GT.1.)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_HIJING(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_HIJING(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   110     KFL1D=KFL1A+INT(RLU_HIJING(0)+0.5)*(KFL1B-KFL1A) 
104           KFL1E=KFL1A+KFL1B-KFL1D   
105           IF((KFL1D.EQ.3.AND.RLU_HIJING(0).GT.PARDM).OR.(KFL1E.EQ.3.AND.   
106      &    RLU_HIJING(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_HIJING(0))   
111           IF((KFL1E.NE.KFL3A.AND.RLU_HIJING(0).GT.(1.+PAR4M)/MAX(2.,1.
112      $         +PAR4M)).OR.(KFL1E.EQ.KFL3A.AND.RLU_HIJING(0).GT.2./MAX(2
113      $         .,1.+PAR4M)))GOTO 110  
114           KFLDS=3   
115           IF(KFL1E.NE.KFL3A) KFLDS=2*INT(RLU_HIJING(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_HIJING(0)) 
125         IF(KFLA.EQ.3) KMUL=INT(PARJ(12)+RLU_HIJING(0)) 
126         IF(KFLA.GE.4) KMUL=INT(PARJ(13)+RLU_HIJING(0)) 
127         IF(KMUL.EQ.0.AND.PARJ(14).GT.0.) THEN   
128           IF(RLU_HIJING(0).LT.PARJ(14)) KMUL=2 
129         ELSEIF(KMUL.EQ.1.AND.PARJ(15)+PARJ(16)+PARJ(17).GT.0.) THEN 
130           RMUL=RLU_HIJING(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_HIJING(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...Generate diquark flavour.   
151       ELSE  
152   120   IF(KF1A.LE.10.AND.KF2A.EQ.0) THEN   
153           KFLA=KF1A 
154   130     KFLB=1+INT((2.+PAR2*PAR3)*RLU_HIJING(0)) 
155           KFLC=1+INT((2.+PAR2*PAR3)*RLU_HIJING(0)) 
156           KFLDS=1   
157           IF(KFLB.GE.KFLC) KFLDS=3  
158           IF(KFLDS.EQ.1.AND.PAR4*RLU_HIJING(0).GT.1.) GOTO 130 
159           IF(KFLDS.EQ.3.AND.PAR4.LT.RLU_HIJING(0)) GOTO 130    
160           KFL3=ISIGN(1000*MAX(KFLB,KFLC)+100*MIN(KFLB,KFLC)+KFLDS,KFL1) 
161     
162 C...Take diquark flavour from input.    
163         ELSEIF(KF1A.LE.10) THEN 
164           KFLA=KF1A 
165           KFLB=MOD(KF2A/1000,10)    
166           KFLC=MOD(KF2A/100,10) 
167           KFLDS=MOD(KF2A,10)    
168     
169 C...Generate (or take from input) quark to go with diquark. 
170         ELSE    
171           IF(KF2A.EQ.0) KFL3=ISIGN(1+INT((2.+PAR2)*RLU_HIJING(0)),KFL1)    
172           KFLA=KF2A+IABS(KFL3)  
173           KFLB=MOD(KF1A/1000,10)    
174           KFLC=MOD(KF1A/100,10) 
175           KFLDS=MOD(KF1A,10)    
176         ENDIF   
177     
178 C...SU(6) factors for formation of baryon. Try again if fails.  
179         KBARY=KFLDS 
180         IF(KFLDS.EQ.3.AND.KFLB.NE.KFLC) KBARY=5 
181         IF(KFLA.NE.KFLB.AND.KFLA.NE.KFLC) KBARY=KBARY+1 
182         WT=PARF(60+KBARY)+PARJ(18)*PARF(70+KBARY)   
183         IF(MBARY.EQ.1.AND.MSTJ(12).GE.2) THEN   
184           WTDQ=PARS0    
185           IF(MAX(KFLB,KFLC).EQ.3) WTDQ=PARS1    
186           IF(MIN(KFLB,KFLC).EQ.3) WTDQ=PARS2    
187           IF(KFLDS.EQ.1) WTDQ=WTDQ/(3.*PAR4M)   
188           IF(KFLDS.EQ.1) WT=WT*(1.+WTDQ)/(1.+PARSM/(3.*PAR4M))  
189           IF(KFLDS.EQ.3) WT=WT*(1.+WTDQ)/(1.+PARSM) 
190         ENDIF   
191         IF(KF2A.EQ.0.AND.WT.LT.RLU_HIJING(0)) GOTO 120 
192     
193 C...Form baryon. Distinguish Lambda- and Sigmalike baryons. 
194         KFLD=MAX(KFLA,KFLB,KFLC)    
195         KFLF=MIN(KFLA,KFLB,KFLC)    
196         KFLE=KFLA+KFLB+KFLC-KFLD-KFLF   
197         KFLS=2  
198         IF((PARF(60+KBARY)+PARJ(18)*PARF(70+KBARY))*RLU_HIJING(0).GT.  
199      &  PARF(60+KBARY)) KFLS=4  
200         KFLL=0  
201         IF(KFLS.EQ.2.AND.KFLD.GT.KFLE.AND.KFLE.GT.KFLF) THEN    
202           IF(KFLDS.EQ.1.AND.KFLA.EQ.KFLD) KFLL=1    
203           IF(KFLDS.EQ.1.AND.KFLA.NE.KFLD) KFLL=INT(0.25+RLU_HIJING(0)) 
204           IF(KFLDS.EQ.3.AND.KFLA.NE.KFLD) KFLL=INT(0.75+RLU_HIJING(0)) 
205         ENDIF   
206         IF(KFLL.EQ.0) KF=ISIGN(1000*KFLD+100*KFLE+10*KFLF+KFLS,KFL1)    
207         IF(KFLL.EQ.1) KF=ISIGN(1000*KFLD+100*KFLF+10*KFLE+KFLS,KFL1)    
208       ENDIF 
209       RETURN    
210     
211 C...Use tabulated probabilities to select new flavour and hadron.   
212   140 IF(KTAB2.EQ.0.AND.MSTJ(12).LE.0) THEN 
213         KT3L=1  
214         KT3U=6  
215       ELSEIF(KTAB2.EQ.0.AND.KTAB1.GE.7.AND.MSTJ(12).LE.1) THEN  
216         KT3L=1  
217         KT3U=6  
218       ELSEIF(KTAB2.EQ.0) THEN   
219         KT3L=1  
220         KT3U=22 
221       ELSE  
222         KT3L=KTAB2  
223         KT3U=KTAB2  
224       ENDIF 
225       RFL=0.    
226       DO 150 KTS=0,2    
227       DO 150 KT3=KT3L,KT3U  
228       RFL=RFL+PARF(120+80*KTAB1+25*KTS+KT3) 
229   150 CONTINUE  
230       RFL=RLU_HIJING(0)*RFL    
231       DO 160 KTS=0,2    
232       KTABS=KTS 
233       DO 160 KT3=KT3L,KT3U  
234       KTAB3=KT3 
235       RFL=RFL-PARF(120+80*KTAB1+25*KTS+KT3) 
236   160 IF(RFL.LE.0.) GOTO 170    
237   170 CONTINUE  
238     
239 C...Reconstruct flavour of produced quark/diquark.  
240       IF(KTAB3.LE.6) THEN   
241         KFL3A=KTAB3 
242         KFL3B=0 
243         KFL3=ISIGN(KFL3A,KFL1*(2*KTAB1-13)) 
244       ELSE  
245         KFL3A=1 
246         IF(KTAB3.GE.8) KFL3A=2  
247         IF(KTAB3.GE.11) KFL3A=3 
248         IF(KTAB3.GE.16) KFL3A=4 
249         KFL3B=(KTAB3-6-KFL3A*(KFL3A-2))/2   
250         KFL3=1000*KFL3A+100*KFL3B+1 
251         IF(KFL3A.EQ.KFL3B.OR.KTAB3.NE.6+KFL3A*(KFL3A-2)+2*KFL3B) KFL3=  
252      &  KFL3+2  
253         KFL3=ISIGN(KFL3,KFL1*(13-2*KTAB1))  
254       ENDIF 
255     
256 C...Reconstruct meson code. 
257       IF(KFL3A.EQ.KFL1A.AND.KFL3B.EQ.KFL1B.AND.(KFL3A.LE.3.OR.  
258      &KFL3B.NE.0)) THEN 
259          RFL=RLU_HIJING(0)*(PARF(143+80*KTAB1+25*KTABS)+PARF(144+80
260      $        *KTAB1+25*KTABS)+PARF(145+80*KTAB1+25*KTABS))  
261         KF=110+2*KTABS+1    
262         IF(RFL.GT.PARF(143+80*KTAB1+25*KTABS)) KF=220+2*KTABS+1 
263         IF(RFL.GT.PARF(143+80*KTAB1+25*KTABS)+PARF(144+80*KTAB1+    
264      &  25*KTABS)) KF=330+2*KTABS+1 
265       ELSEIF(KTAB1.LE.6.AND.KTAB3.LE.6) THEN    
266         KFLA=MAX(KTAB1,KTAB3)   
267         KFLB=MIN(KTAB1,KTAB3)   
268         KFS=ISIGN(1,KFL1)   
269         IF(KFLA.NE.KF1A) KFS=-KFS   
270         KF=(100*KFLA+10*KFLB+2*KTABS+1)*KFS*(-1)**KFLA  
271       ELSEIF(KTAB1.GE.7.AND.KTAB3.GE.7) THEN    
272         KFS=ISIGN(1,KFL1)   
273         IF(KFL1A.EQ.KFL3A) THEN 
274           KFLA=MAX(KFL1B,KFL3B) 
275           KFLB=MIN(KFL1B,KFL3B) 
276           IF(KFLA.NE.KFL1B) KFS=-KFS    
277         ELSEIF(KFL1A.EQ.KFL3B) THEN 
278           KFLA=KFL3A    
279           KFLB=KFL1B    
280           KFS=-KFS  
281         ELSEIF(KFL1B.EQ.KFL3A) THEN 
282           KFLA=KFL1A    
283           KFLB=KFL3B    
284         ELSEIF(KFL1B.EQ.KFL3B) THEN 
285           KFLA=MAX(KFL1A,KFL3A) 
286           KFLB=MIN(KFL1A,KFL3A) 
287           IF(KFLA.NE.KFL1A) KFS=-KFS    
288         ELSE    
289            CALL LUERRM_HIJING(2
290      $          ,'(LUKFDI_HIJING:) 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_HIJING(KF) 
316       IF(KC.EQ.0) THEN  
317          CALL LUERRM_HIJING(2
318      $        ,'(LUKFDI_HIJING:) user-defined flavour probabilities'
319      $        //'failed')   
320         GOTO 100    
321       ENDIF 
322     
323       RETURN    
324       END