pyhidiff.f

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

    
C*********************************************************************  
    
      SUBROUTINE pyhidiff
    
C...Handles diffractive and elastic scattering. 
      common/lujets/n,k(9000,5),p(9000,5),v(9000,5)
      SAVE /lujets/ 
      common/ludat1/mstu(200),paru(200),mstj(200),parj(200) 
      SAVE /ludat1/ 
      common/pyhipars/mstp(200),parp(200),msti(200),pari(200) 
      SAVE /pyhipars/ 
      common/pyhiint1/mint(400),vint(400) 
      SAVE /pyhiint1/ 
    
C...Reset K, P and V vectors. Store incoming particles. 
      DO 100 jt=1,mstp(126)+10  
      i=mint(83)+jt 
      DO 100 j=1,5  
      k(i,j)=0  
      p(i,j)=0. 
  100 v(i,j)=0. 
      n=mint(84)    
      mint(3)=0 
      mint(21)=0    
      mint(22)=0    
      mint(23)=0    
      mint(24)=0    
      mint(4)=4 
      DO 110 jt=1,2 
      i=mint(83)+jt 
      k(i,1)=21 
      k(i,2)=mint(10+jt)    
      p(i,5)=vint(2+jt) 
      p(i,3)=vint(5)*(-1)**(jt+1)   
  110 p(i,4)=sqrt(p(i,3)**2+p(i,5)**2)  
      mint(6)=2 
    
C...Subprocess; kinematics. 
      isub=mint(1)  
      sqlam=(vint(2)-vint(63)-vint(64))**2-4.*vint(63)*vint(64) 
      pz=sqrt(sqlam)/(2.*vint(1))   
      DO 150 jt=1,2 
      i=mint(83)+jt 
      pe=(vint(2)+vint(62+jt)-vint(65-jt))/(2.*vint(1)) 
    
C...Elastically scattered particle. 
      IF(mint(16+jt).LE.0) THEN 
        n=n+1   
        k(n,1)=1    
        k(n,2)=k(i,2)   
        k(n,3)=i+2  
        p(n,3)=pz*(-1)**(jt+1)  
        p(n,4)=pe   
        p(n,5)=p(i,5)   
    
C...Diffracted particle: valence quark kicked out.  
      ELSEIF(mstp(101).EQ.1) THEN   
        n=n+2   
        k(n-1,1)=2  
        k(n,1)=1    
        k(n-1,3)=i+2    
        k(n,3)=i+2  
        CALL pyhispli(k(i,2),21,k(n,2),k(n-1,2))  
        p(n-1,5)=ulmass(k(n-1,2))   
        p(n,5)=ulmass(k(n,2))   
        sqlam=(vint(62+jt)-p(n-1,5)**2-p(n,5)**2)**2-   
     &  4.*p(n-1,5)**2*p(n,5)**2    
        p(n-1,3)=(pe*sqrt(sqlam)+pz*(vint(62+jt)+p(n-1,5)**2-   
     &  p(n,5)**2))/(2.*vint(62+jt))*(-1)**(jt+1)   
        p(n-1,4)=sqrt(p(n-1,3)**2+p(n-1,5)**2)  
        p(n,3)=pz*(-1)**(jt+1)-p(n-1,3) 
        p(n,4)=sqrt(p(n,3)**2+p(n,5)**2)    
    
C...Diffracted particle: gluon kicked out.  
      ELSE  
        n=n+3   
        k(n-2,1)=2  
        k(n-1,1)=2  
        k(n,1)=1    
        k(n-2,3)=i+2    
        k(n-1,3)=i+2    
        k(n,3)=i+2  
        CALL pyhispli(k(i,2),21,k(n,2),k(n-2,2))  
        k(n-1,2)=21 
        p(n-2,5)=ulmass(k(n-2,2))   
        p(n-1,5)=0. 
        p(n,5)=ulmass(k(n,2))   
C...Energy distribution for particle into two jets. 
  120   imb=1   
        IF(mod(k(i,2)/1000,10).NE.0) imb=2  
        chik=parp(92+2*imb) 
        IF(mstp(92).LE.1) THEN  
          IF(imb.EQ.1) chi=rlu(0)   
          IF(imb.EQ.2) chi=1.-sqrt(rlu(0))  
        ELSEIF(mstp(92).EQ.2) THEN  
          chi=1.-rlu(0)**(1./(1.+chik)) 
        ELSEIF(mstp(92).EQ.3) THEN  
          cut=2.*0.3/vint(1)    
  130     chi=rlu(0)**2 
          IF((chi**2/(chi**2+cut**2))**0.25*(1.-chi)**chik.LT.  
     &    rlu(0)) goto 130  
        ELSE    
          cut=2.*0.3/vint(1)    
          cutr=(1.+sqrt(1.+cut**2))/cut 
  140     chir=cut*cutr**rlu(0) 
          chi=(chir**2-cut**2)/(2.*chir)    
          IF((1.-chi)**chik.LT.rlu(0)) goto 140 
        ENDIF   
        IF(chi.LT.p(n,5)**2/vint(62+jt).OR.chi.GT.1.-p(n-2,5)**2/   
     &  vint(62+jt)) goto 120   
        sqm=p(n-2,5)**2/(1.-chi)+p(n,5)**2/chi  
        IF((sqrt(sqm)+parj(32))**2.GE.vint(62+jt)) goto 120 
        pzi=(pe*(vint(62+jt)-sqm)+pz*(vint(62+jt)+sqm))/    
     &  (2.*vint(62+jt))    
        pei=sqrt(pzi**2+sqm)    
        pqqp=(1.-chi)*(pei+pzi) 
        p(n-2,3)=0.5*(pqqp-p(n-2,5)**2/pqqp)*(-1)**(jt+1)   
        p(n-2,4)=sqrt(p(n-2,3)**2+p(n-2,5)**2)  
        p(n-1,3)=(pz-pzi)*(-1)**(jt+1)  
        p(n-1,4)=abs(p(n-1,3))  
        p(n,3)=pzi*(-1)**(jt+1)-p(n-2,3)    
        p(n,4)=sqrt(p(n,3)**2+p(n,5)**2)    
      ENDIF 
    
C...Documentation lines.    
      k(i+2,1)=21   
      IF(mint(16+jt).EQ.0) k(i+2,2)=mint(10+jt) 
      IF(mint(16+jt).NE.0) k(i+2,2)=10*(mint(10+jt)/10) 
      k(i+2,3)=i    
      p(i+2,3)=pz*(-1)**(jt+1)  
      p(i+2,4)=pe   
      p(i+2,5)=sqrt(vint(62+jt))    
  150 CONTINUE  
    
C...Rotate outgoing partons/particles using cos(theta). 
      CALL ludbrb(mint(83)+3,n,acos(vint(23)),vint(24),0d0,0d0,0d0) 
    
      RETURN    
      END