More volume overlaps corrected
[u/mrichter/AliRoot.git] / ISAJET / code / heavyx.F
CommitLineData
0795afa3 1#include "isajet/pilot.h"
2 SUBROUTINE HEAVYX(X,EPS)
3C
4C GENERATE X FOR HEAVY PARTICLE FRAGMENTATION ACCORDING TO
5C THE PETERSON FORM
6C D(X)=1/(X*(1-1/X-EPS/(1-X))**2)
7C =D0(X)*D1(X)*D2(X)
8C D0(X)=(1-X)**2/((1-X)**2+EPS)**2
9C D1(X)=X
10C D2(X)=(((1-X)**2+EPS)/((1-X)**2+EPS*X))**2
11C USING X=1-Y**POW
12C
13 DATA ALN4/1.3863/
14C
15C CHOOSE POW FOR X=1-Y**POW.
16C 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
29C
30C GENERATE Z ACCORDING TO (1-Z)**2/((1-Z)**2+EPS*Z)**2
311 CONTINUE
32 Y=RANF()
33 Z=1.-Y**POW
34C
35 D0Z=(1.-Z)**2/((1.-Z)**2+EPS)**2*POW*Y**(POW-1.)
36 IF(D0Z.LT.RANF()*D0MX) GO TO 1
37C
38C 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
42C
43C GOOD X
44 X=Z
45 RETURN
46 END