pyxkfl.f

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

 
C*********************************************************************
 
C...PYXKFL
C...Selects flavour for produced qqbar pair.
 
      SUBROUTINE pyxkfl(KFL,ECM,ECMC,KFLC)
 
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...Calculate maximum weight in QED or QFD case.
      IF(mstj(102).LE.1) THEN
        rfmax=4d0/9d0
      ELSE
        poll=1d0-parj(131)*parj(132)
        sff=1d0/(16d0*paru(102)*(1d0-paru(102)))
        sfw=ecmc**4/((ecmc**2-parj(123)**2)**2+(parj(123)*parj(124))**2)
        sfi=sfw*(1d0-(parj(123)/ecmc)**2)
        ve=4d0*paru(102)-1d0
        hf1i=sfi*sff*(ve*poll+parj(132)-parj(131))
        hf1w=sfw*sff**2*((ve**2+1d0)*poll+2d0*ve*(parj(132)-parj(131)))
        rfmax=max(4d0/9d0*poll-4d0/3d0*(1d0-8d0*paru(102)/3d0)*hf1i+
     &  ((1d0-8d0*paru(102)/3d0)**2+1d0)*hf1w,1d0/9d0*poll+2d0/3d0*
     &  (-1d0+4d0*paru(102)/3d0)*hf1i+((-1d0+4d0*paru(102)/3d0)**2+
     &  1d0)*hf1w)
      ENDIF
 
C...Choose flavour. Gives charge and velocity.
      ntry=0
  100 ntry=ntry+1
      IF(ntry.GT.100) THEN
        CALL pyerrm(14,'(PYXKFL:) caught in an infinite loop')
        kflc=0
        RETURN
      ENDIF
      kflc=kfl
      IF(kfl.LE.0) kflc=1+int(mstj(104)*pyr(0))
      mstj(93)=1
      pmq=pymass(kflc)
      IF(ecm.LT.2d0*pmq+parj(127)) goto 100
      qf=kchg(kflc,1)/3d0
      vq=1d0
      IF(mod(mstj(103),2).EQ.1) vq=sqrt(max(0d0,1d0-(2d0*pmq/ecmc)**2))
 
C...Calculate weight in QED or QFD case.
      IF(mstj(102).LE.1) THEN
        rf=qf**2
        rfv=0.5d0*vq*(3d0-vq**2)*qf**2
      ELSE
        vf=sign(1d0,qf)-4d0*qf*paru(102)
        rf=qf**2*poll-2d0*qf*vf*hf1i+(vf**2+1d0)*hf1w
        rfv=0.5d0*vq*(3d0-vq**2)*(qf**2*poll-2d0*qf*vf*hf1i+vf**2*hf1w)+
     &  vq**3*hf1w
        IF(rfv.GT.0d0) parj(171)=min(1d0,vq**3*hf1w/rfv)
      ENDIF
 
C...Weighting or new event (radiative photon). Cross-section update.
      IF(kfl.LE.0.AND.rf.LT.pyr(0)*rfmax) goto 100
      parj(158)=parj(158)+1d0
      IF(ecmc.LT.2d0*pmq+parj(127).OR.rfv.LT.pyr(0)*rf) kflc=0
      IF(mstj(107).LE.0.AND.kflc.EQ.0) goto 100
      IF(kflc.NE.0) parj(159)=parj(159)+1d0
      parj(144)=parj(157)*parj(159)/parj(158)
      parj(148)=parj(144)*86.8d0/ecm**2
 
      RETURN
      END