New functions and constructors added and some other fixes.
[u/mrichter/AliRoot.git] / PYTHIA / jetset / luzdis.F
1  
2 C********************************************************************* 
3  
4       SUBROUTINE LUZDIS(KFL1,KFL2,PR,Z) 
5  
6 C...Purpose: to generate the longitudinal splitting variable z. 
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...Check if heavy flavour fragmentation. 
12       KFLA=IABS(KFL1) 
13       KFLB=IABS(KFL2) 
14       KFLH=KFLA 
15       IF(KFLA.GE.10) KFLH=MOD(KFLA/1000,10) 
16  
17 C...Lund symmetric scaling function: determine parameters of shape. 
18       IF(MSTJ(11).EQ.1.OR.(MSTJ(11).EQ.3.AND.KFLH.LE.3).OR. 
19      &MSTJ(11).GE.4) THEN 
20         FA=PARJ(41) 
21         IF(MSTJ(91).EQ.1) FA=PARJ(43) 
22         IF(KFLB.GE.10) FA=FA+PARJ(45) 
23         FBB=PARJ(42) 
24         IF(MSTJ(91).EQ.1) FBB=PARJ(44) 
25         FB=FBB*PR 
26         FC=1. 
27         IF(KFLA.GE.10) FC=FC-PARJ(45) 
28         IF(KFLB.GE.10) FC=FC+PARJ(45) 
29         IF(MSTJ(11).GE.4.AND.KFLH.GE.4.AND.KFLH.LE.5) THEN 
30           FRED=PARJ(46) 
31           IF(MSTJ(11).EQ.5.AND.KFLH.EQ.5) FRED=PARJ(47) 
32           FC=FC+FRED*FBB*PARF(100+KFLH)**2 
33         ELSEIF(MSTJ(11).GE.4.AND.KFLH.GE.6.AND.KFLH.LE.8) THEN 
34           FRED=PARJ(46) 
35           IF(MSTJ(11).EQ.5) FRED=PARJ(48) 
36           FC=FC+FRED*FBB*PMAS(KFLH,1)**2 
37         ENDIF 
38         MC=1 
39         IF(ABS(FC-1.).GT.0.01) MC=2 
40  
41 C...Determine position of maximum. Special cases for a = 0 or a = c. 
42         IF(FA.LT.0.02) THEN 
43           MA=1 
44           ZMAX=1. 
45           IF(FC.GT.FB) ZMAX=FB/FC 
46         ELSEIF(ABS(FC-FA).LT.0.01) THEN 
47           MA=2 
48           ZMAX=FB/(FB+FC) 
49         ELSE 
50           MA=3 
51           ZMAX=0.5*(FB+FC-SQRT((FB-FC)**2+4.*FA*FB))/(FC-FA) 
52           IF(ZMAX.GT.0.9999.AND.FB.GT.100.) ZMAX=MIN(ZMAX,1.-FA/FB) 
53         ENDIF 
54  
55 C...Subdivide z range if distribution very peaked near endpoint. 
56         MMAX=2 
57         IF(ZMAX.LT.0.1) THEN 
58           MMAX=1 
59           ZDIV=2.75*ZMAX 
60           IF(MC.EQ.1) THEN 
61             FINT=1.-LOG(ZDIV) 
62           ELSE 
63             ZDIVC=ZDIV**(1.-FC) 
64             FINT=1.+(1.-1./ZDIVC)/(FC-1.) 
65           ENDIF 
66         ELSEIF(ZMAX.GT.0.85.AND.FB.GT.1.) THEN 
67           MMAX=3 
68           FSCB=SQRT(4.+(FC/FB)**2) 
69           ZDIV=FSCB-1./ZMAX-(FC/FB)*LOG(ZMAX*0.5*(FSCB+FC/FB)) 
70           IF(MA.GE.2) ZDIV=ZDIV+(FA/FB)*LOG(1.-ZMAX) 
71           ZDIV=MIN(ZMAX,MAX(0.,ZDIV)) 
72           FINT=1.+FB*(1.-ZDIV) 
73         ENDIF 
74  
75 C...Choice of z, preweighted for peaks at low or high z. 
76   100   Z=RLU(0) 
77         FPRE=1. 
78         IF(MMAX.EQ.1) THEN 
79           IF(FINT*RLU(0).LE.1.) THEN 
80             Z=ZDIV*Z 
81           ELSEIF(MC.EQ.1) THEN 
82             Z=ZDIV**Z 
83             FPRE=ZDIV/Z 
84           ELSE 
85             Z=1./(ZDIVC+Z*(1.-ZDIVC))**(1./(1.-FC)) 
86             FPRE=(ZDIV/Z)**FC 
87           ENDIF 
88         ELSEIF(MMAX.EQ.3) THEN 
89           IF(FINT*RLU(0).LE.1.) THEN 
90             Z=ZDIV+LOG(Z)/FB 
91             FPRE=EXP(FB*(Z-ZDIV)) 
92           ELSE 
93             Z=ZDIV+Z*(1.-ZDIV) 
94           ENDIF 
95         ENDIF 
96  
97 C...Weighting according to correct formula. 
98         IF(Z.LE.0..OR.Z.GE.1.) GOTO 100 
99         FEXP=FC*LOG(ZMAX/Z)+FB*(1./ZMAX-1./Z) 
100         IF(MA.GE.2) FEXP=FEXP+FA*LOG((1.-Z)/(1.-ZMAX)) 
101         FVAL=EXP(MAX(-50.,MIN(50.,FEXP))) 
102         IF(FVAL.LT.RLU(0)*FPRE) GOTO 100 
103  
104 C...Generate z according to Field-Feynman, SLAC, (1-z)**c OR z**c. 
105       ELSE 
106         FC=PARJ(50+MAX(1,KFLH)) 
107         IF(MSTJ(91).EQ.1) FC=PARJ(59) 
108   110   Z=RLU(0) 
109         IF(FC.GE.0..AND.FC.LE.1.) THEN 
110           IF(FC.GT.RLU(0)) Z=1.-Z**(1./3.) 
111         ELSEIF(FC.GT.-1.AND.FC.LT.0.) THEN 
112           IF(-4.*FC*Z*(1.-Z)**2.LT.RLU(0)*((1.-Z)**2-FC*Z)**2) GOTO 110 
113         ELSE 
114           IF(FC.GT.0.) Z=1.-Z**(1./FC) 
115           IF(FC.LT.0.) Z=Z**(-1./FC) 
116         ENDIF 
117       ENDIF 
118  
119       RETURN 
120       END