luzdis.f

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

    
C*********************************************************************  
    
      SUBROUTINE luzdis(KFL1,KFL2,PR,Z) 
    
C...Purpose: to generate the longitudinal splitting variable z. 
      common/ludat1/mstu(200),paru(200),mstj(200),parj(200) 
      SAVE /ludat1/ 
    
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)) THEN   
        fa=parj(41) 
        IF(mstj(91).EQ.1) fa=parj(43)   
        IF(kflb.GE.10) fa=fa+parj(45)   
        fb=parj(42)*pr  
        IF(mstj(91).EQ.1) fb=parj(44)*pr    
        fc=1.   
        IF(kfla.GE.10) fc=fc-parj(45)   
        IF(kflb.GE.10) fc=fc+parj(45)   
        mc=1    
        IF(abs(fc-1.).GT.0.01) mc=2 
    
C...Determine position of maximum. Special cases for a = 0 or a = c.    
        IF(fa.LT.0.02) THEN 
          ma=1  
          zmax=1.   
          IF(fc.GT.fb) zmax=fb/fc   
        ELSEIF(abs(fc-fa).LT.0.01) THEN 
          ma=2  
          zmax=fb/(fb+fc)   
        ELSE    
          ma=3  
          zmax=0.5*(fb+fc-sqrt((fb-fc)**2+4.*fa*fb))/(fc-fa)    
          IF(zmax.GT.0.99.AND.fb.GT.100.) zmax=1.-fa/fb 
        ENDIF   
    
C...Subdivide z range if distribution very peaked near endpoint.    
        mmax=2  
        IF(zmax.LT.0.1) THEN    
          mmax=1    
          zdiv=2.75*zmax    
          IF(mc.EQ.1) THEN  
            fint=1.-log(zdiv)   
          ELSE  
            zdivc=zdiv**(1.-fc) 
            fint=1.+(1.-1./zdivc)/(fc-1.)   
          ENDIF 
        ELSEIF(zmax.GT.0.85.AND.fb.GT.1.) THEN  
          mmax=3    
          fscb=sqrt(4.+(fc/fb)**2)  
          zdiv=fscb-1./zmax-(fc/fb)*log(zmax*0.5*(fscb+fc/fb))  
          IF(ma.GE.2) zdiv=zdiv+(fa/fb)*log(1.-zmax)    
          zdiv=min(zmax,max(0.,zdiv))   
          fint=1.+fb*(1.-zdiv)  
        ENDIF   
    
C...Choice of z, preweighted for peaks at low or high z.    
  100   z=rlu(0)    
        fpre=1. 
        IF(mmax.EQ.1) THEN  
          IF(fint*rlu(0).LE.1.) THEN    
            z=zdiv*z    
          ELSEIF(mc.EQ.1) THEN  
            z=zdiv**z   
            fpre=zdiv/z 
          ELSE  
            z=1./(zdivc+z*(1.-zdivc))**(1./(1.-fc)) 
            fpre=(zdiv/z)**fc   
          ENDIF 
        ELSEIF(mmax.EQ.3) THEN  
          IF(fint*rlu(0).LE.1.) THEN    
            z=zdiv+log(z)/fb    
            fpre=exp(fb*(z-zdiv))   
          ELSE  
            z=zdiv+z*(1.-zdiv)  
          ENDIF 
        ENDIF   
    
C...Weighting according to correct formula. 
        IF(z.LE.fb/(50.+fb).OR.z.GE.1.) goto 100    
        fval=(zmax/z)**fc*exp(fb*(1./zmax-1./z))    
        IF(ma.GE.2) fval=((1.-z)/(1.-zmax))**fa*fval    
        IF(fval.LT.rlu(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=rlu(0)    
        IF(fc.GE.0..AND.fc.LE.1.) THEN  
          IF(fc.GT.rlu(0)) z=1.-z**(1./3.)  
        ELSEIF(fc.GT.-1.) THEN  
          IF(-4.*fc*z*(1.-z)**2.LT.rlu(0)*((1.-z)**2-fc*z)**2) goto 110 
        ELSE    
          IF(fc.GT.0.) z=1.-z**(1./fc)  
          IF(fc.LT.0.) z=z**(-1./fc)    
        ENDIF   
      ENDIF 
    
      RETURN    
      END