luxkfl.f

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

    
C*********************************************************************  
    
      SUBROUTINE luxkfl(KFL,ECM,ECMC,KFLC)  
    
C...Purpose: to select flavour for produced qqbar pair. 
      common/ludat1/mstu(200),paru(200),mstj(200),parj(200) 
      SAVE /ludat1/ 
      common/ludat2/kchg(500,3),pmas(500,4),parf(2000),vckm(4,4)    
      SAVE /ludat2/ 
    
C...Calculate maximum weight in QED or QFD case.    
      IF(mstj(102).LE.1) THEN   
        rfmax=4./9. 
      ELSE  
        poll=1.-parj(131)*parj(132) 
        sff=1./(16.*paru(102)*(1.-paru(102)))   
        sfw=ecmc**4/((ecmc**2-parj(123)**2)**2+(parj(123)*parj(124))**2)    
        sfi=sfw*(1.-(parj(123)/ecmc)**2)    
        ve=4.*paru(102)-1.  
        hf1i=sfi*sff*(ve*poll+parj(132)-parj(131))  
        hf1w=sfw*sff**2*((ve**2+1.)*poll+2.*ve*(parj(132)-parj(131)))   
        rfmax=max(4./9.*poll-4./3.*(1.-8.*paru(102)/3.)*hf1i+   
     &  ((1.-8.*paru(102)/3.)**2+1.)*hf1w,1./9.*poll+2./3.* 
     &  (-1.+4.*paru(102)/3.)*hf1i+((-1.+4.*paru(102)/3.)**2+1.)*hf1w)  
      ENDIF 
    
C...Choose flavour. Gives charge and velocity.  
      ntry=0    
  100 ntry=ntry+1   
      IF(ntry.GT.100) THEN  
        CALL luerrm(14,'(LUXKFL:) caught in an infinite loop')  
        kflc=0  
        RETURN  
      ENDIF 
      kflc=kfl  
      IF(kfl.LE.0) kflc=1+int(mstj(104)*rlu(0)) 
      mstj(93)=1    
      pmq=ulmass(kflc)  
      IF(ecm.LT.2.*pmq+parj(127)) goto 100  
      qf=kchg(kflc,1)/3.    
      vq=1. 
      IF(mod(mstj(103),2).EQ.1) vq=sqrt(max(0.,1.-(2.*pmq/ecmc)**2))    
    
C...Calculate weight in QED or QFD case.    
      IF(mstj(102).LE.1) THEN   
        rf=qf**2    
        rfv=0.5*vq*(3.-vq**2)*qf**2 
      ELSE  
        vf=sign(1.,qf)-4.*qf*paru(102)  
        rf=qf**2*poll-2.*qf*vf*hf1i+(vf**2+1.)*hf1w 
        rfv=0.5*vq*(3.-vq**2)*(qf**2*poll-2.*qf*vf*hf1i+vf**2*hf1w)+    
     &  vq**3*hf1w  
      ENDIF 
    
C...Weighting or new event (radiative photon). Cross-section update.    
      IF(kfl.LE.0.AND.rf.LT.rlu(0)*rfmax) goto 100  
      parj(158)=parj(158)+1.    
      IF(ecmc.LT.2.*pmq+parj(127).OR.rfv.LT.rlu(0)*rf) kflc=0   
      IF(mstj(107).LE.0.AND.kflc.EQ.0) goto 100 
      IF(kflc.NE.0) parj(159)=parj(159)+1.  
      parj(144)=parj(157)*parj(159)/parj(158)   
      parj(148)=parj(144)*86.8/ecm**2   
    
      RETURN    
      END