pythiaeRHIC/pythia/pyrobo.f

0001
0002 C*********************************************************************
0003
0004 C...PYROBO
0005 C...Performs rotations and boosts.
0006
0007       SUBROUTINE PYROBO(IMI,IMA,THE,PHI,BEX,BEY,BEZ)
0008
0009 C...Double precision and integer declarations.
0010       IMPLICIT DOUBLE PRECISION(A-H, O-Z)
0011       IMPLICIT INTEGER(I-N)
0012       INTEGER PYK,PYCHGE,PYCOMP
0013 C...Commonblocks.
0014       COMMON/PYJETS/N,NPAD,K(4000,5),P(4000,5),V(4000,5)
0015       COMMON/PYDAT1/MSTU(200),PARU(200),MSTJ(200),PARJ(200)
0016       SAVE /PYJETS/,/PYDAT1/
0017 C...Local arrays.
0018       DIMENSION ROT(3,3),PR(3),VR(3),DP(4),DV(4)
0019
0020 C...Find and check range of rotation/boost.
0021       IMIN=IMI
0022       IF(IMIN.LE.0) IMIN=1
0023       IF(MSTU(1).GT.0) IMIN=MSTU(1)
0024       IMAX=IMA
0025       IF(IMAX.LE.0) IMAX=N
0026       IF(MSTU(2).GT.0) IMAX=MSTU(2)
0027       IF(IMIN.GT.MSTU(4).OR.IMAX.GT.MSTU(4)) THEN
0028         CALL PYERRM(11,'(PYROBO:) range outside PYJETS memory')
0029         RETURN
0030       ENDIF
0031
0032 C...Optional resetting of V (when not set before.)
0033       IF(MSTU(33).NE.0) THEN
0034         DO 110 I=MIN(IMIN,MSTU(4)),MIN(IMAX,MSTU(4))
0035           DO 100 J=1,5
0036             V(I,J)=0D0
0037   100     CONTINUE
0038   110   CONTINUE
0039         MSTU(33)=0
0040       ENDIF
0041
0042 C...Rotate, typically from z axis to direction (theta,phi).
0043       IF(THE**2+PHI**2.GT.1D-20) THEN
0044         ROT(1,1)=COS(THE)*COS(PHI)
0045         ROT(1,2)=-SIN(PHI)
0046         ROT(1,3)=SIN(THE)*COS(PHI)
0047         ROT(2,1)=COS(THE)*SIN(PHI)
0048         ROT(2,2)=COS(PHI)
0049         ROT(2,3)=SIN(THE)*SIN(PHI)
0050         ROT(3,1)=-SIN(THE)
0051         ROT(3,2)=0D0
0052         ROT(3,3)=COS(THE)
0053         DO 140 I=IMIN,IMAX
0054           IF(K(I,1).LE.0) GOTO 140
0055           DO 120 J=1,3
0056             PR(J)=P(I,J)
0057             VR(J)=V(I,J)
0058   120     CONTINUE
0059           DO 130 J=1,3
0060             P(I,J)=ROT(J,1)*PR(1)+ROT(J,2)*PR(2)+ROT(J,3)*PR(3)
0061             V(I,J)=ROT(J,1)*VR(1)+ROT(J,2)*VR(2)+ROT(J,3)*VR(3)
0062   130     CONTINUE
0063   140   CONTINUE
0064       ENDIF
0065
0066 C...Boost, typically from rest to momentum/energy=beta.
0067       IF(BEX**2+BEY**2+BEZ**2.GT.1D-20) THEN
0068         DBX=BEX
0069         DBY=BEY
0070         DBZ=BEZ
0071         DB=SQRT(DBX**2+DBY**2+DBZ**2)
0072         EPS1=1D0-1D-12
0073         IF(DB.GT.EPS1) THEN
0074 C...Rescale boost vector if too close to unity.
0075           CALL PYERRM(3,'(PYROBO:) boost vector too large')
0076           DBX=DBX*(EPS1/DB)
0077           DBY=DBY*(EPS1/DB)
0078           DBZ=DBZ*(EPS1/DB)
0079           DB=EPS1
0080         ENDIF
0081         DGA=1D0/SQRT(1D0-DB**2)
0082         DO 160 I=IMIN,IMAX
0083           IF(K(I,1).LE.0) GOTO 160
0084           DO 150 J=1,4
0085             DP(J)=P(I,J)
0086             DV(J)=V(I,J)
0087   150     CONTINUE
0088           DBP=DBX*DP(1)+DBY*DP(2)+DBZ*DP(3)
0089           DGABP=DGA*(DGA*DBP/(1D0+DGA)+DP(4))
0090           P(I,1)=DP(1)+DGABP*DBX
0091           P(I,2)=DP(2)+DGABP*DBY
0092           P(I,3)=DP(3)+DGABP*DBZ
0093           P(I,4)=DGA*(DP(4)+DBP)
0094           DBV=DBX*DV(1)+DBY*DV(2)+DBZ*DV(3)
0095           DGABV=DGA*(DGA*DBV/(1D0+DGA)+DV(4))
0096           V(I,1)=DV(1)+DGABV*DBX
0097           V(I,2)=DV(2)+DGABV*DBY
0098           V(I,3)=DV(3)+DGABV*DBZ
0099           V(I,4)=DGA*(DV(4)+DBV)
0100   160   CONTINUE
0101       ENDIF
0102
0103       RETURN
0104       END