pythiaeRHIC/pythia/pyjurf.f

0001
0002 C*********************************************************************
0003
0004 C...PYJURF
0005 C...From three given input vectors in PJU the boost VJU from
0006 C...the "lab frame" to the junction rest frame is constructed.
0007
0008       SUBROUTINE PYJURF(PJU,VJU)
0009
0010 C...Double precision and integer declarations.
0011       IMPLICIT DOUBLE PRECISION(A-H, O-Z)
0012       IMPLICIT INTEGER(I-N)
0013
0014 C...Input, output and local arrays.
0015       DIMENSION PJU(3,5),VJU(5),PSUM(5),A(3,3),PENEW(3),PCM(5,5)
0016       DATA TWOPI/6.283186D0/
0017
0018 C...Calculate masses and other invariants.
0019       DO 100 J=1,4
0020         PSUM(J)=PJU(1,J)+PJU(2,J)+PJU(3,J)
0021   100 CONTINUE
0022       PSUM2=PSUM(4)**2-PSUM(1)**2-PSUM(2)**2-PSUM(3)**2
0023       PSUM(5)=SQRT(PSUM2)
0024       DO 120 I=1,3
0025         DO 110 J=1,3
0026           A(I,J)=PJU(I,4)*PJU(J,4)-PJU(I,1)*PJU(J,1)-
0027      &    PJU(I,2)*PJU(J,2)-PJU(I,3)*PJU(J,3)
0028   110   CONTINUE
0029   120 CONTINUE
0030
0031 C...Pick I to be most massive parton and J to be the one closest to I.
0032       ITRY=0
0033       I=1
0034       IF(A(2,2).GT.A(1,1)) I=2
0035       IF(A(3,3).GT.MAX(A(1,1),A(2,2))) I=3
0036   130 ITRY=ITRY+1
0037       J=1+MOD(I,3)
0038       K=1+MOD(J,3)
0039       IF(A(I,K)**2*A(J,J).LT.A(I,J)**2*A(K,K)) THEN
0040         K=1+MOD(I,3)
0041         J=1+MOD(K,3)
0042       ENDIF
0043       PMI2=A(I,I)
0044       PMJ2=A(J,J)
0045       PMK2=A(K,K)
0046       AIJ=A(I,J)
0047       AIK=A(I,K)
0048       AJK=A(J,K)
0049
0050 C...Trivial find new parton energies if all three partons are massless.
0051       IF(PMI2.LT.1D-4) THEN
0052         PEI=SQRT(2D0*AIK*AIJ/(3D0*AJK))
0053         PEJ=SQRT(2D0*AJK*AIJ/(3D0*AIK))
0054         PEK=SQRT(2D0*AIK*AJK/(3D0*AIJ))
0055
0056 C...Else find momentum range for parton I and values at extremes.
0057       ELSE
0058         PAIMIN=0D0
0059         PEIMIN=SQRT(PMI2)
0060         PEJMIN=AIJ/PEIMIN
0061         PEKMIN=AIK/PEIMIN
0062         PAJMIN=SQRT(MAX(0D0,PEJMIN**2-PMJ2))
0063         PAKMIN=SQRT(MAX(0D0,PEKMIN**2-PMK2))
0064         FMIN=PEJMIN*PEKMIN+0.5D0*PAJMIN*PAKMIN-AJK
0065         PEIMAX=(AIJ+AIK)/SQRT(PMJ2+PMK2+2D0*AJK)
0066         IF(PMJ2.GT.1D-4) PEIMAX=AIJ/SQRT(PMJ2)
0067         PAIMAX=SQRT(MAX(0D0,PEIMAX**2-PMI2))
0068         HI=PEIMAX**2-0.25D0*PAIMAX**2
0069         PAJMAX=(PEIMAX*SQRT(MAX(0D0,AIJ**2-PMJ2*HI))-
0070      &  0.5D0*PAIMAX*AIJ)/HI
0071         PAKMAX=(PEIMAX*SQRT(MAX(0D0,AIK**2-PMK2*HI))-
0072      &  0.5D0*PAIMAX*AIK)/HI
0073         PEJMAX=SQRT(PAJMAX**2+PMJ2)
0074         PEKMAX=SQRT(PAKMAX**2+PMK2)
0075         FMAX=PEJMAX*PEKMAX+0.5D0*PAJMAX*PAKMAX-AJK
0076
0077 C...If unexpected values at upper endpoint then pick another parton.
0078         IF(FMAX.GT.0D0.AND.ITRY.LE.2) THEN
0079           I1=1+MOD(I,3)
0080           IF(A(I1,I1).GE.1D-4) THEN
0081             I=I1
0082             GOTO 130
0083           ENDIF
0084           ITRY=ITRY+1
0085           I1=1+MOD(I,3)
0086           IF(ITRY.LE.2.AND.A(I1,I1).GE.1D-4) THEN
0087             I=I1
0088             GOTO 130
0089           ENDIF
0090         ENDIF
0091
0092 C..Start binary + linear search to find solution inside range.
0093         ITER=0
0094         ITMIN=0
0095         ITMAX=0
0096         PAI=0.5D0*(PAIMIN+PAIMAX)
0097   140   ITER=ITER+1
0098
0099 C...Derive momentum of other two partons and distance to root.
0100         PEI=SQRT(PAI**2+PMI2)
0101         HI=PEI**2-0.25D0*PAI**2
0102         PAJ=(PEI*SQRT(MAX(0D0,AIJ**2-PMJ2*HI))-0.5D0*PAI*AIJ)/HI
0103         PEJ=SQRT(PAJ**2+PMJ2)
0104         PAK=(PEI*SQRT(MAX(0D0,AIK**2-PMK2*HI))-0.5D0*PAI*AIK)/HI
0105         PEK=SQRT(PAK**2+PMK2)
0106         FNOW=PEJ*PEK+0.5D0*PAJ*PAK-AJK
0107
0108 C...Pick next I momentum to explore, hopefully closer to root.
0109         IF(FNOW.GT.0D0) THEN
0110           PAIMIN=PAI
0111           FMIN=FNOW
0112           ITMIN=ITMIN+1
0113         ELSE
0114           PAIMAX=PAI
0115           FMAX=FNOW
0116           ITMAX=ITMAX+1
0117         ENDIF
0118         IF((ITER.LT.10.OR.ITMIN.LE.1.OR.ITMAX.LE.1).AND.ITER.LT.20)
0119      &  THEN
0120           PAI=0.5D0*(PAIMIN+PAIMAX)
0121           GOTO 140
0122         ELSEIF(ITER.LT.40.AND.FMIN.GT.0D0.AND.FMAX.LT.0D0.AND.
0123      &  ABS(FNOW).GT.1D-12*PSUM2) THEN
0124           PAI=PAIMIN+(PAIMAX-PAIMIN)*FMIN/(FMIN-FMAX)
0125           GOTO 140
0126         ENDIF
0127       ENDIF
0128
0129 C...Now know energies in junction rest frame.
0130       PENEW(I)=PEI
0131       PENEW(J)=PEJ
0132       PENEW(K)=PEK
0133
0134 C...Boost (copy of) partons to their rest frame.
0135       VXCM=-PSUM(1)/PSUM(5)
0136       VYCM=-PSUM(2)/PSUM(5)
0137       VZCM=-PSUM(3)/PSUM(5)
0138       GAMCM=SQRT(1D0+VXCM**2+VYCM**2+VZCM**2)
0139       DO 150 I=1,3
0140         FAC1=PJU(I,1)*VXCM+PJU(I,2)*VYCM+PJU(I,3)*VZCM
0141         FAC2=FAC1/(1D0+GAMCM)+PJU(I,4)
0142         PCM(I,1)=PJU(I,1)+FAC2*VXCM
0143         PCM(I,2)=PJU(I,2)+FAC2*VYCM
0144         PCM(I,3)=PJU(I,3)+FAC2*VZCM
0145         PCM(I,4)=PJU(I,4)*GAMCM+FAC1
0146         PCM(I,5)=SQRT(PCM(I,1)**2+PCM(I,2)**2+PCM(I,3)**2)
0147   150 CONTINUE
0148
0149 C...Construct difference vectors and boost to junction rest frame.
0150       DO 160 J=1,3
0151         PCM(4,J)=PCM(1,J)/PCM(1,4)-PCM(2,J)/PCM(2,4)
0152         PCM(5,J)=PCM(1,J)/PCM(1,4)-PCM(3,J)/PCM(3,4)
0153   160 CONTINUE
0154       PCM(4,4)=PENEW(1)/PCM(1,4)-PENEW(2)/PCM(2,4)
0155       PCM(5,4)=PENEW(1)/PCM(1,4)-PENEW(3)/PCM(3,4)
0156       PCM4S=PCM(4,1)**2+PCM(4,2)**2+PCM(4,3)**2
0157       PCM5S=PCM(5,1)**2+PCM(5,2)**2+PCM(5,3)**2
0158       PCM45=PCM(4,1)*PCM(5,1)+PCM(4,2)*PCM(5,2)+PCM(4,3)*PCM(5,3)
0159       C4=(PCM5S*PCM(4,4)-PCM45*PCM(5,4))/(PCM4S*PCM5S-PCM45**2)
0160       C5=(PCM4S*PCM(5,4)-PCM45*PCM(4,4))/(PCM4S*PCM5S-PCM45**2)
0161       VXJU=C4*PCM(4,1)+C5*PCM(5,1)
0162       VYJU=C4*PCM(4,2)+C5*PCM(5,2)
0163       VZJU=C4*PCM(4,3)+C5*PCM(5,3)
0164       GAMJU=SQRT(1D0+VXJU**2+VYJU**2+VZJU**2)
0165
0166 C...Add two boosts, giving final result.
0167       FCM=(VXJU*VXCM+VYJU*VYCM+VZJU*VZCM)/(1+GAMCM)+GAMJU
0168       VJU(1)=VXJU+FCM*VXCM
0169       VJU(2)=VYJU+FCM*VYCM
0170       VJU(3)=VZJU+FCM*VZCM
0171       VJU(4)=SQRT(1D0+VJU(1)**2+VJU(2)**2+VJU(3)**2)
0172       VJU(5)=1D0
0173
0174 C...In case of error in reconstruction: revert to CM frame of system.
0175       CTH12=(PCM(1,1)*PCM(2,1)+PCM(1,2)*PCM(2,2)+PCM(1,3)*PCM(2,3))/
0176      &(PCM(1,5)*PCM(2,5))
0177       CTH13=(PCM(1,1)*PCM(3,1)+PCM(1,2)*PCM(3,2)+PCM(1,3)*PCM(3,3))/
0178      &(PCM(1,5)*PCM(3,5))
0179       CTH23=(PCM(2,1)*PCM(3,1)+PCM(2,2)*PCM(3,2)+PCM(2,3)*PCM(3,3))/
0180      &(PCM(2,5)*PCM(3,5))
0181       ERRCCM=(CTH12+0.5D0)**2+(CTH13+0.5D0)**2+(CTH23+0.5D0)**2
0182       ERRTCM=TWOPI-ACOS(CTH12)-ACOS(CTH13)-ACOS(CTH23)
0183       DO 170 I=1,3
0184         FAC1=PJU(I,1)*VJU(1)+PJU(I,2)*VJU(2)+PJU(I,3)*VJU(3)
0185         FAC2=FAC1/(1D0+VJU(4))+PJU(I,4)
0186         PCM(I,1)=PJU(I,1)+FAC2*VJU(1)
0187         PCM(I,2)=PJU(I,2)+FAC2*VJU(2)
0188         PCM(I,3)=PJU(I,3)+FAC2*VJU(3)
0189         PCM(I,4)=PJU(I,4)*VJU(4)+FAC1
0190         PCM(I,5)=SQRT(PCM(I,1)**2+PCM(I,2)**2+PCM(I,3)**2)
0191   170 CONTINUE
0192       CTH12=(PCM(1,1)*PCM(2,1)+PCM(1,2)*PCM(2,2)+PCM(1,3)*PCM(2,3))/
0193      &(PCM(1,5)*PCM(2,5))
0194       CTH13=(PCM(1,1)*PCM(3,1)+PCM(1,2)*PCM(3,2)+PCM(1,3)*PCM(3,3))/
0195      &(PCM(1,5)*PCM(3,5))
0196       CTH23=(PCM(2,1)*PCM(3,1)+PCM(2,2)*PCM(3,2)+PCM(2,3)*PCM(3,3))/
0197      &(PCM(2,5)*PCM(3,5))
0198       ERRCJU=(CTH12+0.5D0)**2+(CTH13+0.5D0)**2+(CTH23+0.5D0)**2
0199       ERRTJU=TWOPI-ACOS(CTH12)-ACOS(CTH13)-ACOS(CTH23)
0200       IF(ERRCJU+ERRTJU.GT.ERRCCM+ERRTCM) THEN
0201         VJU(1)=VXCM
0202         VJU(2)=VYCM
0203         VJU(3)=VZCM
0204         VJU(4)=GAMCM
0205       ENDIF
0206
0207       RETURN
0208       END