pyrobo.f

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

 
C*********************************************************************
 
C...PYROBO
C...Performs rotations and boosts.
 
      SUBROUTINE pyrobo(IMI,IMA,THE,PHI,BEX,BEY,BEZ)
 
C...Double precision and integer declarations.
      IMPLICIT DOUBLE PRECISION(a-h, o-z)
      IMPLICIT INTEGER(i-n)
      INTEGER pyk,pychge,pycomp
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.
      dimension rot(3,3),pr(3),vr(3),dp(4),dv(4)
 
C...Find and check range of rotation/boost.
      imin=imi
      IF(imin.LE.0) imin=1
      IF(mstu(1).GT.0) imin=mstu(1)
      imax=ima
      IF(imax.LE.0) imax=n
      IF(mstu(2).GT.0) imax=mstu(2)
      IF(imin.GT.mstu(4).OR.imax.GT.mstu(4)) THEN
        CALL pyerrm(11,'(PYROBO:) range outside PYJETS memory')
        RETURN
      ENDIF
 
C...Optional resetting of V (when not set before.)
      IF(mstu(33).NE.0) THEN
        DO 110 i=min(imin,mstu(4)),min(imax,mstu(4))
          DO 100 j=1,5
            v(i,j)=0d0
  100     CONTINUE
  110   CONTINUE
        mstu(33)=0
      ENDIF
 
C...Rotate, typically from z axis to direction (theta,phi).
      IF(the**2+phi**2.GT.1d-20) THEN
        rot(1,1)=cos(the)*cos(phi)
        rot(1,2)=-sin(phi)
        rot(1,3)=sin(the)*cos(phi)
        rot(2,1)=cos(the)*sin(phi)
        rot(2,2)=cos(phi)
        rot(2,3)=sin(the)*sin(phi)
        rot(3,1)=-sin(the)
        rot(3,2)=0d0
        rot(3,3)=cos(the)
        DO 140 i=imin,imax
          IF(k(i,1).LE.0) goto 140
          DO 120 j=1,3
            pr(j)=p(i,j)
            vr(j)=v(i,j)
  120     CONTINUE
          DO 130 j=1,3
            p(i,j)=rot(j,1)*pr(1)+rot(j,2)*pr(2)+rot(j,3)*pr(3)
            v(i,j)=rot(j,1)*vr(1)+rot(j,2)*vr(2)+rot(j,3)*vr(3)
  130     CONTINUE
  140   CONTINUE
      ENDIF
 
C...Boost, typically from rest to momentum/energy=beta.
      IF(bex**2+bey**2+bez**2.GT.1d-20) THEN
        dbx=bex
        dby=bey
        dbz=bez
        db=sqrt(dbx**2+dby**2+dbz**2)
        eps1=1d0-1d-12
        IF(db.GT.eps1) THEN
C...Rescale boost vector if too close to unity.
          CALL pyerrm(3,'(PYROBO:) boost vector too large')
          dbx=dbx*(eps1/db)
          dby=dby*(eps1/db)
          dbz=dbz*(eps1/db)
          db=eps1
        ENDIF
        dga=1d0/sqrt(1d0-db**2)
        DO 160 i=imin,imax
          IF(k(i,1).LE.0) goto 160
          DO 150 j=1,4
            dp(j)=p(i,j)
            dv(j)=v(i,j)
  150     CONTINUE
          dbp=dbx*dp(1)+dby*dp(2)+dbz*dp(3)
          dgabp=dga*(dga*dbp/(1d0+dga)+dp(4))
          p(i,1)=dp(1)+dgabp*dbx
          p(i,2)=dp(2)+dgabp*dby
          p(i,3)=dp(3)+dgabp*dbz
          p(i,4)=dga*(dp(4)+dbp)
          dbv=dbx*dv(1)+dby*dv(2)+dbz*dv(3)
          dgabv=dga*(dga*dbv/(1d0+dga)+dv(4))
          v(i,1)=dv(1)+dgabv*dbx
          v(i,2)=dv(2)+dgabv*dby
          v(i,3)=dv(3)+dgabv*dbz
          v(i,4)=dga*(dv(4)+dbv)
  160   CONTINUE
      ENDIF
 
      RETURN
      END