luradk.f

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

    
C*********************************************************************  
    
      SUBROUTINE luradk(ECM,MK,PAK,THEK,PHIK,ALPK)  
    
C...Purpose: to generate initial state photon radiation.    
      common/ludat1/mstu(200),paru(200),mstj(200),parj(200) 
      SAVE /ludat1/ 
    
C...Function: cumulative hard photon spectrum in QFD case.  
      fxk(xx)=2.*log(xx)+parj(161)*log(1.-xx)+parj(162)*xx+ 
     &parj(163)*log((xx-szm)**2+szw**2)+parj(164)*atan((xx-szm)/szw)    
    
C...Determine whether radiative photon or not.  
      mk=0  
      pak=0.    
      IF(parj(160).LT.rlu(0)) RETURN    
      mk=1  
    
C...Photon energy range. Find photon momentum in QED case.  
      xkl=parj(135) 
      xku=min(parj(136),1.-(2.*parj(127)/ecm)**2)   
      IF(mstj(102).LE.1) THEN   
  100   xk=1./(1.+(1./xkl-1.)*((1./xku-1.)/(1./xkl-1.))**rlu(0))    
        IF(1.+(1.-xk)**2.LT.2.*rlu(0)) goto 100 
    
C...Ditto in QFD case, by numerical inversion of integrated spectrum.   
      ELSE  
        szm=1.-(parj(123)/ecm)**2   
        szw=parj(123)*parj(124)/ecm**2  
        fxkl=fxk(xkl)   
        fxku=fxk(xku)   
        fxkd=1e-4*(fxku-fxkl)   
        fxkr=fxkl+rlu(0)*(fxku-fxkl)    
        nxk=0   
  110   nxk=nxk+1   
        xk=0.5*(xkl+xku)    
        fxkv=fxk(xk)    
        IF(fxkv.GT.fxkr) THEN   
          xku=xk    
          fxku=fxkv 
        ELSE    
          xkl=xk    
          fxkl=fxkv 
        ENDIF   
        IF(nxk.LT.15.AND.fxku-fxkl.GT.fxkd) goto 110    
        xk=xkl+(xku-xkl)*(fxkr-fxkl)/(fxku-fxkl)    
      ENDIF 
      pak=0.5*ecm*xk    
    
C...Photon polar and azimuthal angle.   
      pme=2.*(ulmass(11)/ecm)**2    
  120 cthm=pme*(2./pme)**rlu(0) 
      IF(1.-(xk**2*cthm*(1.-0.5*cthm)+2.*(1.-xk)*pme/max(pme,   
     &cthm*(1.-0.5*cthm)))/(1.+(1.-xk)**2).LT.rlu(0)) goto 120  
      cthe=1.-cthm  
      IF(rlu(0).GT.0.5) cthe=-cthe  
      sthe=sqrt(max(0.,(cthm-pme)*(2.-cthm)))   
      thek=ulangl(cthe,sthe)    
      phik=paru(2)*rlu(0)   
    
C...Rotation angle for hadronic system. 
      sgn=1.    
      IF(0.5*(2.-xk*(1.-cthe))**2/((2.-xk)**2+(xk*cthe)**2).GT. 
     &rlu(0)) sgn=-1.   
      alpk=asin(sgn*sthe*(xk-sgn*(2.*sqrt(1.-xk)-2.+xk)*cthe)/  
     &(2.-xk*(1.-sgn*cthe)))    
    
      RETURN    
      END