pyzdis.f

Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file pyzdis.f

 
C*********************************************************************
 
C...PYZDIS
C...Generates the longitudinal splitting variable z.
 
      SUBROUTINE pyzdis(KFL1,KFL2,PR,Z)
 
C...Double precision and integer declarations.
      IMPLICIT DOUBLE PRECISION(a-h, o-z)
      IMPLICIT INTEGER(i-n)
      INTEGER pyk,pychge,pycomp
C...Commonblocks.
      common/pydat1/mstu(200),paru(200),mstj(200),parj(200)
      common/pydat2/kchg(500,4),pmas(500,4),parf(2000),vckm(4,4)
      SAVE /pydat1/,/pydat2/
 
C...Check if heavy flavour fragmentation.
      kfla=iabs(kfl1)
      kflb=iabs(kfl2)
      kflh=kfla
      IF(kfla.GE.10) kflh=mod(kfla/1000,10)
 
C...Lund symmetric scaling function: determine parameters of shape.
      IF(mstj(11).EQ.1.OR.(mstj(11).EQ.3.AND.kflh.LE.3).OR.
     &mstj(11).GE.4) THEN
        fa=parj(41)
        IF(mstj(91).EQ.1) fa=parj(43)
        IF(kflb.GE.10) fa=fa+parj(45)
        fbb=parj(42)
        IF(mstj(91).EQ.1) fbb=parj(44)
        fb=fbb*pr
        fc=1d0
        IF(kfla.GE.10) fc=fc-parj(45)
        IF(kflb.GE.10) fc=fc+parj(45)
        IF(mstj(11).GE.4.AND.(kflh.EQ.4.OR.kflh.EQ.5)) THEN
          fred=parj(46)
          IF(mstj(11).EQ.5.AND.kflh.EQ.5) fred=parj(47)
          fc=fc+fred*fbb*parf(100+kflh)**2
        ENDIF
        mc=1
        IF(abs(fc-1d0).GT.0.01d0) mc=2
 
C...Determine position of maximum. Special cases for a = 0 or a = c.
        IF(fa.LT.0.02d0) THEN
          ma=1
          zmax=1d0
          IF(fc.GT.fb) zmax=fb/fc
        ELSEIF(abs(fc-fa).LT.0.01d0) THEN
          ma=2
          zmax=fb/(fb+fc)
        ELSE
          ma=3
          zmax=0.5d0*(fb+fc-sqrt((fb-fc)**2+4d0*fa*fb))/(fc-fa)
          IF(zmax.GT.0.9999d0.AND.fb.GT.100d0) zmax=min(zmax,1d0-fa/fb)
        ENDIF
 
C...Subdivide z range if distribution very peaked near endpoint.
        mmax=2
        IF(zmax.LT.0.1d0) THEN
          mmax=1
          zdiv=2.75d0*zmax
          IF(mc.EQ.1) THEN
            fint=1d0-log(zdiv)
          ELSE
            zdivc=zdiv**(1d0-fc)
            fint=1d0+(1d0-1d0/zdivc)/(fc-1d0)
          ENDIF
        ELSEIF(zmax.GT.0.85d0.AND.fb.GT.1d0) THEN
          mmax=3
          fscb=sqrt(4d0+(fc/fb)**2)
          zdiv=fscb-1d0/zmax-(fc/fb)*log(zmax*0.5d0*(fscb+fc/fb))
          IF(ma.GE.2) zdiv=zdiv+(fa/fb)*log(1d0-zmax)
          zdiv=min(zmax,max(0d0,zdiv))
          fint=1d0+fb*(1d0-zdiv)
        ENDIF
 
C...Choice of z, preweighted for peaks at low or high z.
  100   z=pyr(0)
        fpre=1d0
        IF(mmax.EQ.1) THEN
          IF(fint*pyr(0).LE.1d0) THEN
            z=zdiv*z
          ELSEIF(mc.EQ.1) THEN
            z=zdiv**z
            fpre=zdiv/z
          ELSE
            z=(zdivc+z*(1d0-zdivc))**(1d0/(1d0-fc))
            fpre=(zdiv/z)**fc
          ENDIF
        ELSEIF(mmax.EQ.3) THEN
          IF(fint*pyr(0).LE.1d0) THEN
            z=zdiv+log(z)/fb
            fpre=exp(fb*(z-zdiv))
          ELSE
            z=zdiv+z*(1d0-zdiv)
          ENDIF
        ENDIF
 
C...Weighting according to correct formula.
        IF(z.LE.0d0.OR.z.GE.1d0) goto 100
        fexp=fc*log(zmax/z)+fb*(1d0/zmax-1d0/z)
        IF(ma.GE.2) fexp=fexp+fa*log((1d0-z)/(1d0-zmax))
        fval=exp(max(-50d0,min(50d0,fexp)))
        IF(fval.LT.pyr(0)*fpre) goto 100
 
C...Generate z according to Field-Feynman, SLAC, (1-z)**c OR z**c.
      ELSE
        fc=parj(50+max(1,kflh))
        IF(mstj(91).EQ.1) fc=parj(59)
  110   z=pyr(0)
        IF(fc.GE.0d0.AND.fc.LE.1d0) THEN
          IF(fc.GT.pyr(0)) z=1d0-z**(1d0/3d0)
        ELSEIF(fc.GT.-1.AND.fc.LT.0d0) THEN
          IF(-4d0*fc*z*(1d0-z)**2.LT.pyr(0)*((1d0-z)**2-fc*z)**2)
     &    goto 110
        ELSE
          IF(fc.GT.0d0) z=1d0-z**(1d0/fc)
          IF(fc.LT.0d0) z=z**(-1d0/fc)
        ENDIF
      ENDIF
 
      RETURN
      END