EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
luzdis.f
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file 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  SAVE /ludat1/
9 
10 C...Check if heavy flavour fragmentation.
11  kfla=iabs(kfl1)
12  kflb=iabs(kfl2)
13  kflh=kfla
14  IF(kfla.GE.10) kflh=mod(kfla/1000,10)
15 
16 C...Lund symmetric scaling function: determine parameters of shape.
17  IF(mstj(11).EQ.1.OR.(mstj(11).EQ.3.AND.kflh.LE.3)) THEN
18  fa=parj(41)
19  IF(mstj(91).EQ.1) fa=parj(43)
20  IF(kflb.GE.10) fa=fa+parj(45)
21  fb=parj(42)*pr
22  IF(mstj(91).EQ.1) fb=parj(44)*pr
23  fc=1.
24  IF(kfla.GE.10) fc=fc-parj(45)
25  IF(kflb.GE.10) fc=fc+parj(45)
26  mc=1
27  IF(abs(fc-1.).GT.0.01) mc=2
28 
29 C...Determine position of maximum. Special cases for a = 0 or a = c.
30  IF(fa.LT.0.02) THEN
31  ma=1
32  zmax=1.
33  IF(fc.GT.fb) zmax=fb/fc
34  ELSEIF(abs(fc-fa).LT.0.01) THEN
35  ma=2
36  zmax=fb/(fb+fc)
37  ELSE
38  ma=3
39  zmax=0.5*(fb+fc-sqrt((fb-fc)**2+4.*fa*fb))/(fc-fa)
40  IF(zmax.GT.0.99.AND.fb.GT.100.) zmax=1.-fa/fb
41  ENDIF
42 
43 C...Subdivide z range if distribution very peaked near endpoint.
44  mmax=2
45  IF(zmax.LT.0.1) THEN
46  mmax=1
47  zdiv=2.75*zmax
48  IF(mc.EQ.1) THEN
49  fint=1.-log(zdiv)
50  ELSE
51  zdivc=zdiv**(1.-fc)
52  fint=1.+(1.-1./zdivc)/(fc-1.)
53  ENDIF
54  ELSEIF(zmax.GT.0.85.AND.fb.GT.1.) THEN
55  mmax=3
56  fscb=sqrt(4.+(fc/fb)**2)
57  zdiv=fscb-1./zmax-(fc/fb)*log(zmax*0.5*(fscb+fc/fb))
58  IF(ma.GE.2) zdiv=zdiv+(fa/fb)*log(1.-zmax)
59  zdiv=min(zmax,max(0.,zdiv))
60  fint=1.+fb*(1.-zdiv)
61  ENDIF
62 
63 C...Choice of z, preweighted for peaks at low or high z.
64  100 z=rlu(0)
65  fpre=1.
66  IF(mmax.EQ.1) THEN
67  IF(fint*rlu(0).LE.1.) THEN
68  z=zdiv*z
69  ELSEIF(mc.EQ.1) THEN
70  z=zdiv**z
71  fpre=zdiv/z
72  ELSE
73  z=1./(zdivc+z*(1.-zdivc))**(1./(1.-fc))
74  fpre=(zdiv/z)**fc
75  ENDIF
76  ELSEIF(mmax.EQ.3) THEN
77  IF(fint*rlu(0).LE.1.) THEN
78  z=zdiv+log(z)/fb
79  fpre=exp(fb*(z-zdiv))
80  ELSE
81  z=zdiv+z*(1.-zdiv)
82  ENDIF
83  ENDIF
84 
85 C...Weighting according to correct formula.
86  IF(z.LE.fb/(50.+fb).OR.z.GE.1.) goto 100
87  fval=(zmax/z)**fc*exp(fb*(1./zmax-1./z))
88  IF(ma.GE.2) fval=((1.-z)/(1.-zmax))**fa*fval
89  IF(fval.LT.rlu(0)*fpre) goto 100
90 
91 C...Generate z according to Field-Feynman, SLAC, (1-z)**c OR z**c.
92  ELSE
93  fc=parj(50+max(1,kflh))
94  IF(mstj(91).EQ.1) fc=parj(59)
95  110 z=rlu(0)
96  IF(fc.GE.0..AND.fc.LE.1.) THEN
97  IF(fc.GT.rlu(0)) z=1.-z**(1./3.)
98  ELSEIF(fc.GT.-1.) THEN
99  IF(-4.*fc*z*(1.-z)**2.LT.rlu(0)*((1.-z)**2-fc*z)**2) goto 110
100  ELSE
101  IF(fc.GT.0.) z=1.-z**(1./fc)
102  IF(fc.LT.0.) z=z**(-1./fc)
103  ENDIF
104  ENDIF
105 
106  RETURN
107  END