hirobo.f

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

    
C*********************************************************************  
C THIS SUBROUTINE IS ONLY FOR THE USE OF HIJING TO ROTATE OR BOOST
C       THE FOUR MOMENTUM ONLY
C*********************************************************************
    
      SUBROUTINE hirobo(THE,PHI,BEX,BEY,BEZ)    
    
C...Purpose: to perform rotations and boosts.   
      IMPLICIT DOUBLE PRECISION(d)  
      common/lujets/n,k(9000,5),p(9000,5),v(9000,5)
      SAVE /lujets/ 
      common/ludat1/mstu(200),paru(200),mstj(200),parj(200) 
      SAVE /ludat1/ 
      dimension rot(3,3),pr(3),vr(3),dp(4),dv(4)    
    
C...Find range of rotation/boost. Convert boost to double precision.    
      imin=1    
      IF(mstu(1).GT.0) imin=mstu(1) 
      imax=n    
      IF(mstu(2).GT.0) imax=mstu(2) 
      dbx=bex   
      dby=bey   
      dbz=bez   
    
C...Check range of rotation/boost.  
      IF(imin.GT.mstu(4).OR.imax.GT.mstu(4)) THEN   
        CALL luerrm(11,'(LUROBO:) range outside LUJETS memory') 
        RETURN  
      ENDIF 
    
C...Rotate, typically from z axis to direction (theta,phi). 
      IF(the**2+phi**2.GT.1e-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)=0. 
        rot(3,3)=cos(the)   
        DO 130 i=imin,imax  
        IF(k(i,1).LE.0) goto 130    
        DO 110 j=1,3    
  110   pr(j)=p(i,j)   
        DO 120 j=1,3    
  120   p(i,j)=rot(j,1)*pr(1)+rot(j,2)*pr(2)+rot(j,3)*pr(3) 
  130   CONTINUE    
      ENDIF 
    
C...Boost, typically from rest to momentum/energy=beta. 
      IF(dbx**2+dby**2+dbz**2.GT.1e-20) THEN    
        db=sqrt(dbx**2+dby**2+dbz**2)   
        IF(db.GT.0.99999999d0) THEN 
C...Rescale boost vector if too close to unity. 
          CALL luerrm(3,'(LUROBO:) boost vector too large') 
          dbx=dbx*(0.99999999d0/db) 
          dby=dby*(0.99999999d0/db) 
          dbz=dbz*(0.99999999d0/db) 
          db=0.99999999d0   
        ENDIF   
        dga=1d0/sqrt(1d0-db**2) 
        DO 150 i=imin,imax  
        IF(k(i,1).LE.0) goto 150    
        DO 140 j=1,4    
  140   dp(j)=p(i,j)    
        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)  
  150   CONTINUE    
      ENDIF 
    
      RETURN    
      END