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