pyjurf.f

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

 
C*********************************************************************
 
C...PYJURF
C...From three given input vectors in PJU the boost VJU from
C...the "lab frame" to the junction rest frame is constructed.
 
      SUBROUTINE pyjurf(PJU,VJU)
 
C...Double precision and integer declarations.
      IMPLICIT DOUBLE PRECISION(a-h, o-z)
      IMPLICIT INTEGER(i-n)
 
C...Input, output and local arrays.
      dimension pju(3,5),vju(5),psum(5),a(3,3),penew(3),pcm(5,5)
      DATA twopi/6.283186d0/
 
C...Calculate masses and other invariants.
      DO 100 j=1,4
        psum(j)=pju(1,j)+pju(2,j)+pju(3,j)
  100 CONTINUE
      psum2=psum(4)**2-psum(1)**2-psum(2)**2-psum(3)**2
      psum(5)=sqrt(psum2)
      DO 120 i=1,3
        DO 110 j=1,3
          a(i,j)=pju(i,4)*pju(j,4)-pju(i,1)*pju(j,1)-
     &    pju(i,2)*pju(j,2)-pju(i,3)*pju(j,3)
  110   CONTINUE
  120 CONTINUE
 
C...Pick I to be most massive parton and J to be the one closest to I.
      itry=0
      i=1
      IF(a(2,2).GT.a(1,1)) i=2
      IF(a(3,3).GT.max(a(1,1),a(2,2))) i=3
  130 itry=itry+1
      j=1+mod(i,3)
      k=1+mod(j,3)
      IF(a(i,k)**2*a(j,j).LT.a(i,j)**2*a(k,k)) THEN
        k=1+mod(i,3)
        j=1+mod(k,3)
      ENDIF
      pmi2=a(i,i)
      pmj2=a(j,j)
      pmk2=a(k,k)
      aij=a(i,j)
      aik=a(i,k)
      ajk=a(j,k)
 
C...Trivial find new parton energies if all three partons are massless.
      IF(pmi2.LT.1d-4) THEN
        pei=sqrt(2d0*aik*aij/(3d0*ajk))
        pej=sqrt(2d0*ajk*aij/(3d0*aik))
        pek=sqrt(2d0*aik*ajk/(3d0*aij))
 
C...Else find momentum range for parton I and values at extremes.
      ELSE
        paimin=0d0
        peimin=sqrt(pmi2)
        pejmin=aij/peimin
        pekmin=aik/peimin
        pajmin=sqrt(max(0d0,pejmin**2-pmj2))
        pakmin=sqrt(max(0d0,pekmin**2-pmk2))
        fmin=pejmin*pekmin+0.5d0*pajmin*pakmin-ajk
        peimax=(aij+aik)/sqrt(pmj2+pmk2+2d0*ajk)
        IF(pmj2.GT.1d-4) peimax=aij/sqrt(pmj2)
        paimax=sqrt(max(0d0,peimax**2-pmi2))
        hi=peimax**2-0.25d0*paimax**2
        pajmax=(peimax*sqrt(max(0d0,aij**2-pmj2*hi))-
     &  0.5d0*paimax*aij)/hi
        pakmax=(peimax*sqrt(max(0d0,aik**2-pmk2*hi))-
     &  0.5d0*paimax*aik)/hi
        pejmax=sqrt(pajmax**2+pmj2)
        pekmax=sqrt(pakmax**2+pmk2)
        fmax=pejmax*pekmax+0.5d0*pajmax*pakmax-ajk
 
C...If unexpected values at upper endpoint then pick another parton.
        IF(fmax.GT.0d0.AND.itry.LE.2) THEN
          i1=1+mod(i,3)
          IF(a(i1,i1).GE.1d-4) THEN
            i=i1
            goto 130
          ENDIF
          itry=itry+1
          i1=1+mod(i,3)
          IF(itry.LE.2.AND.a(i1,i1).GE.1d-4) THEN
            i=i1
            goto 130
          ENDIF
        ENDIF
 
C..Start binary + linear search to find solution inside range.
        iter=0
        itmin=0
        itmax=0
        pai=0.5d0*(paimin+paimax)
  140   iter=iter+1
 
C...Derive momentum of other two partons and distance to root.
        pei=sqrt(pai**2+pmi2)
        hi=pei**2-0.25d0*pai**2
        paj=(pei*sqrt(max(0d0,aij**2-pmj2*hi))-0.5d0*pai*aij)/hi
        pej=sqrt(paj**2+pmj2)
        pak=(pei*sqrt(max(0d0,aik**2-pmk2*hi))-0.5d0*pai*aik)/hi
        pek=sqrt(pak**2+pmk2)
        fnow=pej*pek+0.5d0*paj*pak-ajk
 
C...Pick next I momentum to explore, hopefully closer to root.
        IF(fnow.GT.0d0) THEN
          paimin=pai
          fmin=fnow
          itmin=itmin+1
        ELSE
          paimax=pai
          fmax=fnow
          itmax=itmax+1
        ENDIF
        IF((iter.LT.10.OR.itmin.LE.1.OR.itmax.LE.1).AND.iter.LT.20)
     &  THEN
          pai=0.5d0*(paimin+paimax)
          goto 140
        ELSEIF(iter.LT.40.AND.fmin.GT.0d0.AND.fmax.LT.0d0.AND.
     &  abs(fnow).GT.1d-12*psum2) THEN
          pai=paimin+(paimax-paimin)*fmin/(fmin-fmax)
          goto 140
        ENDIF
      ENDIF
 
C...Now know energies in junction rest frame.
      penew(i)=pei
      penew(j)=pej
      penew(k)=pek
 
C...Boost (copy of) partons to their rest frame.
      vxcm=-psum(1)/psum(5)
      vycm=-psum(2)/psum(5)
      vzcm=-psum(3)/psum(5)
      gamcm=sqrt(1d0+vxcm**2+vycm**2+vzcm**2)
      DO 150 i=1,3
        fac1=pju(i,1)*vxcm+pju(i,2)*vycm+pju(i,3)*vzcm
        fac2=fac1/(1d0+gamcm)+pju(i,4)
        pcm(i,1)=pju(i,1)+fac2*vxcm
        pcm(i,2)=pju(i,2)+fac2*vycm
        pcm(i,3)=pju(i,3)+fac2*vzcm
        pcm(i,4)=pju(i,4)*gamcm+fac1
        pcm(i,5)=sqrt(pcm(i,1)**2+pcm(i,2)**2+pcm(i,3)**2)
  150 CONTINUE
 
C...Construct difference vectors and boost to junction rest frame.
      DO 160 j=1,3
        pcm(4,j)=pcm(1,j)/pcm(1,4)-pcm(2,j)/pcm(2,4)
        pcm(5,j)=pcm(1,j)/pcm(1,4)-pcm(3,j)/pcm(3,4)
  160 CONTINUE
      pcm(4,4)=penew(1)/pcm(1,4)-penew(2)/pcm(2,4)
      pcm(5,4)=penew(1)/pcm(1,4)-penew(3)/pcm(3,4)
      pcm4s=pcm(4,1)**2+pcm(4,2)**2+pcm(4,3)**2
      pcm5s=pcm(5,1)**2+pcm(5,2)**2+pcm(5,3)**2
      pcm45=pcm(4,1)*pcm(5,1)+pcm(4,2)*pcm(5,2)+pcm(4,3)*pcm(5,3)
      c4=(pcm5s*pcm(4,4)-pcm45*pcm(5,4))/(pcm4s*pcm5s-pcm45**2)
      c5=(pcm4s*pcm(5,4)-pcm45*pcm(4,4))/(pcm4s*pcm5s-pcm45**2)
      vxju=c4*pcm(4,1)+c5*pcm(5,1)
      vyju=c4*pcm(4,2)+c5*pcm(5,2)
      vzju=c4*pcm(4,3)+c5*pcm(5,3)
      gamju=sqrt(1d0+vxju**2+vyju**2+vzju**2)
 
C...Add two boosts, giving final result.
      fcm=(vxju*vxcm+vyju*vycm+vzju*vzcm)/(1+gamcm)+gamju
      vju(1)=vxju+fcm*vxcm
      vju(2)=vyju+fcm*vycm
      vju(3)=vzju+fcm*vzcm
      vju(4)=sqrt(1d0+vju(1)**2+vju(2)**2+vju(3)**2)
      vju(5)=1d0
 
C...In case of error in reconstruction: revert to CM frame of system.
      cth12=(pcm(1,1)*pcm(2,1)+pcm(1,2)*pcm(2,2)+pcm(1,3)*pcm(2,3))/
     &(pcm(1,5)*pcm(2,5))
      cth13=(pcm(1,1)*pcm(3,1)+pcm(1,2)*pcm(3,2)+pcm(1,3)*pcm(3,3))/
     &(pcm(1,5)*pcm(3,5))
      cth23=(pcm(2,1)*pcm(3,1)+pcm(2,2)*pcm(3,2)+pcm(2,3)*pcm(3,3))/
     &(pcm(2,5)*pcm(3,5))
      errccm=(cth12+0.5d0)**2+(cth13+0.5d0)**2+(cth23+0.5d0)**2
      errtcm=twopi-acos(cth12)-acos(cth13)-acos(cth23)
      DO 170 i=1,3
        fac1=pju(i,1)*vju(1)+pju(i,2)*vju(2)+pju(i,3)*vju(3)
        fac2=fac1/(1d0+vju(4))+pju(i,4)
        pcm(i,1)=pju(i,1)+fac2*vju(1)
        pcm(i,2)=pju(i,2)+fac2*vju(2)
        pcm(i,3)=pju(i,3)+fac2*vju(3)
        pcm(i,4)=pju(i,4)*vju(4)+fac1
        pcm(i,5)=sqrt(pcm(i,1)**2+pcm(i,2)**2+pcm(i,3)**2)
  170 CONTINUE
      cth12=(pcm(1,1)*pcm(2,1)+pcm(1,2)*pcm(2,2)+pcm(1,3)*pcm(2,3))/
     &(pcm(1,5)*pcm(2,5))
      cth13=(pcm(1,1)*pcm(3,1)+pcm(1,2)*pcm(3,2)+pcm(1,3)*pcm(3,3))/
     &(pcm(1,5)*pcm(3,5))
      cth23=(pcm(2,1)*pcm(3,1)+pcm(2,2)*pcm(3,2)+pcm(2,3)*pcm(3,3))/
     &(pcm(2,5)*pcm(3,5))
      errcju=(cth12+0.5d0)**2+(cth13+0.5d0)**2+(cth23+0.5d0)**2
      errtju=twopi-acos(cth12)-acos(cth13)-acos(cth23)
      IF(errcju+errtju.GT.errccm+errtcm) THEN
        vju(1)=vxcm
        vju(2)=vycm
        vju(3)=vzcm
        vju(4)=gamcm
      ENDIF
 
      RETURN
      END