More volume overlaps corrected
[u/mrichter/AliRoot.git] / ISAJET / code / heavyx.F
1 #include "isajet/pilot.h"
2       SUBROUTINE HEAVYX(X,EPS)
3 C
4 C          GENERATE X FOR HEAVY PARTICLE FRAGMENTATION ACCORDING TO
5 C          THE PETERSON FORM
6 C          D(X)=1/(X*(1-1/X-EPS/(1-X))**2)
7 C              =D0(X)*D1(X)*D2(X)
8 C          D0(X)=(1-X)**2/((1-X)**2+EPS)**2
9 C          D1(X)=X
10 C          D2(X)=(((1-X)**2+EPS)/((1-X)**2+EPS*X))**2
11 C          USING X=1-Y**POW
12 C
13       DATA ALN4/1.3863/
14 C
15 C          CHOOSE POW FOR X=1-Y**POW.
16 C          GENERATE FLAT IN X IF EPS>1.
17       IF(EPS.LT.1.) THEN
18         POW=ALOG((3.+EPS)/EPS)/ALN4
19         YMX=(EPS*(3.*POW-1.)/(POW+1.))**(.5/POW)
20         ZMX=1-YMX**POW
21         D0MX=(1-ZMX)**2/((1.-ZMX)**2+EPS)**2*POW*YMX**(POW-1.)
22         D2MX=2./(2.-SQRT(EPS))
23       ELSE
24         POW=1.
25         ZMX=0.
26         D0MX=(1.-ZMX)**2/((1.-ZMX)**2+EPS)**2
27         D2MX=1.+EPS
28       ENDIF
29 C
30 C          GENERATE Z ACCORDING TO (1-Z)**2/((1-Z)**2+EPS*Z)**2
31 1     CONTINUE
32       Y=RANF()
33       Z=1.-Y**POW
34 C
35       D0Z=(1.-Z)**2/((1.-Z)**2+EPS)**2*POW*Y**(POW-1.)
36       IF(D0Z.LT.RANF()*D0MX) GO TO 1
37 C
38 C          CHECK REMAINING FACTORS
39       D1=Z
40       D2=(((1.-Z)**2+EPS)/((1.-Z)**2+EPS*Z))**2
41       IF(D1*D2.LT.RANF()*D2MX) GO TO 1
42 C
43 C          GOOD X
44       X=Z
45       RETURN
46       END