pybesq.f

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

 
C*********************************************************************
 
C...PYBESQ
C...Calculates the momentum shift in a system of two particles assuming
C...the relative momentum squared should be shifted to Q2NEW. NI is the
C...last position occupied in /PYJETS/.
 
      SUBROUTINE pybesq(I1,I2,NI,Q2OLD,Q2NEW)
 
C...Double precision and integer declarations.
      IMPLICIT DOUBLE PRECISION(a-h, o-z)
      IMPLICIT INTEGER(i-n)
      INTEGER pyk,pychge,pycomp
C...Parameter statement to help give large particle numbers.
      parameter(ksusy1=1000000,ksusy2=2000000,ktechn=3000000,
     &kexcit=4000000,kdimen=5000000)
C...Commonblocks.
      common/pyjets/n,npad,k(4000,5),p(4000,5),v(4000,5)
      common/pydat1/mstu(200),paru(200),mstj(200),parj(200)
      SAVE /pyjets/,/pydat1/
C...Local arrays and data.
      dimension dp(5)
      SAVE hc1
 
      IF(mstj(55).EQ.0) THEN
        dq2=q2new-q2old
        dp2=(p(i1,1)-p(i2,1))**2+(p(i1,2)-p(i2,2))**2+
     &  (p(i1,3)-p(i2,3))**2
        dp12=p(i1,1)**2+p(i1,2)**2+p(i1,3)**2
     &  -p(i2,1)**2-p(i2,2)**2-p(i2,3)**2
        se=p(i1,4)+p(i2,4)
        de=p(i1,4)-p(i2,4)
        dq2se=dq2+se**2
        da=se*de*dp12-dp2*dq2se
        db=dp2*dq2se-dp12**2
        ha=(da+sqrt(max(da**2+dq2*(dq2+se**2-de**2)*db,0d0)))/(2d0*db)
        DO 100 j=1,3
          pd=ha*(p(i1,j)-p(i2,j))
          p(ni+1,j)=pd
          p(ni+2,j)=-pd
  100   CONTINUE
        RETURN
      ENDIF
 
      k(ni+1,1)=1
      k(ni+2,1)=1
      DO 110 j=1,5
        p(ni+1,j)=p(i1,j)
        p(ni+2,j)=p(i2,j)
        dp(j)=p(i1,j)+p(i2,j)
  110 CONTINUE
 
C...Boost to cms and rotate first particle to z-axis
      CALL pyrobo(ni+1,ni+2,0.0d0,0.0d0,
     &-dp(1)/dp(4),-dp(2)/dp(4),-dp(3)/dp(4))
      phi=pyangl(p(ni+1,1),p(ni+1,2))
      the=pyangl(p(ni+1,3),sqrt(p(ni+1,1)**2+p(ni+1,2)**2))
      s=q2new+(p(i1,5)+p(i2,5))**2
      pz=0.5d0*sqrt(q2new*(s-(p(i1,5)-p(i2,5))**2)/s)
      p(ni+1,1)=0.0d0
      p(ni+1,2)=0.0d0
      p(ni+1,3)=pz
      p(ni+1,4)=sqrt(pz**2+p(i1,5)**2)
      p(ni+2,1)=0.0d0
      p(ni+2,2)=0.0d0
      p(ni+2,3)=-pz
      p(ni+2,4)=sqrt(pz**2+p(i2,5)**2)
      dp(4)=sqrt(dp(1)**2+dp(2)**2+dp(3)**2+s)
      CALL pyrobo(ni+1,ni+2,the,phi,
     &dp(1)/dp(4),dp(2)/dp(4),dp(3)/dp(4))
 
      DO 120 j=1,3
        p(ni+1,j)=p(ni+1,j)-p(i1,j)
        p(ni+2,j)=p(ni+2,j)-p(i2,j)
  120 CONTINUE
 
      RETURN
      END