cmagneticDef.f File Reference

Go to the source code of this file.

Functions/Subroutines
subroutine	cmagneticdef (aTrack, B, leng, dr, dir)

Function/Subroutine Documentation

cmagneticdef()

subroutine cmagneticdef	(	type(track)	aTrack,
		type(magfield)	B,
		real*8	leng,
		type(coord)	dr,
		type(coord)	dir
	)

Definition at line 5 of file cmagneticDef.f.

References cos, cscalerprod(), cvecprod(), d, and d0.

Referenced by cmagdef().

!     ************************************************  
!    This is exact if B is const.
!   Note that:
!      dr must be added  as  newpos = oldpos +  dr
!      dir is the new direction cos.  These two are differenct
!      from the old version.
!
 
      implicit none
#include  "Ztrack.h"


      type(track)::atrack       !input. a charged ptcl
      type(magfield)::b         !input. magnetic field vector.

!                         ptcl and magfiled coord is in Exyz

      real*8  leng                !input. length  travelled in m
      type(coord)::dr           !output. displacement vector in m
      type(coord)::dir          !output. new direction cos. 
!
      type(coord)::tx
      type(coord)::ty
      type(coord)::tz
      type(coord)::w
      type(coord)::r

      real*8 babs, sint, cost, rgt, sinl, cosl
      real*8   pabs, temp
!
!       get abs(B)
      babs = sqrt( b%x**2 + b%y**2 + b%z**2 )
      if(babs .gt. 0.) then
!            p in GeV
         pabs = sqrt(atrack%p%fm%p(4)**2 - atrack%p%mass**2)
         if(pabs .gt. 0.) then
!              transverse gyroradius in m. B in T, p in GeV
!              For negative charge, rgt < 0. This is o.k
            rgt = 3.3358*pabs/atrack%p%charge/babs

            tz%r(1) = b%x/babs
            tz%r(2) = b%y/babs
            tz%r(3) = b%z/babs

            call cscalerprod(atrack%vec%w, b, cost)
            cost =max( min(cost/babs, 1.0d0), -1.d0)
            sint = sqrt(1.0d0-cost**2)
            if(abs(sint) .gt. 1.d-6) then            
               call cvecprod(atrack%vec%w, tz, tx)
               tx%r(1) = tx%r(1)/sint
               tx%r(2) = tx%r(2)/sint
               tx%r(3) = tx%r(3)/sint
               call cvecprod(tz, tx, ty)
!////////////////
!               call checkuv('Tx', sint, Tx)
!               call checkuv('Ty', sint, Ty)
!               call checkuv('Tz', sint, Tz)
!               call checkorth('TxTy', sint, Tx, Ty)
!               call checkorth('TxTz', sint, Tx, Tz)
!//////////////////////               
            else
               tx%r(1) = 1.d0
               tx%r(2) = 0.
               tx%r(3) = 0.
               ty%r(1) = 0.
               ty%r(2) = 1.d0
               ty%r(3) = 0.
            endif
!            get wz0=wz
            call cscalerprod(tz, atrack%vec%w, w%r(3))
            temp = leng /rgt
            sinl = sin(temp)         ! this may be < 0
            cosl = cos(temp)

            w%r(1) = sint*sinl
            w%r(2) = sint*cosl
!///////////
!            call checkuv('w ', sint, w)
!///////////
            r%r(1) = - rgt*sint*(cosl - 1.d0)


            r%r(2) = rgt*sint*sinl
            r%r(3) = leng*w%r(3)
            
!                convert to original system.
            dr%r(1) =
     *          tx%r(1)*r%r(1) + ty%r(1)*r%r(2) + tz%r(1)*r%r(3)
            dr%r(2) =
     *          tx%r(2)*r%r(1) + ty%r(2)*r%r(2) + tz%r(2)*r%r(3)
            dr%r(3) =
     *          tx%r(3)*r%r(1) + ty%r(3)*r%r(2) + tz%r(3)*r%r(3)


            dir%r(1) =
     *          tx%r(1)*w%r(1) + ty%r(1)*w%r(2) + tz%r(1)*w%r(3)
            dir%r(2) =
     *          tx%r(2)*w%r(1) + ty%r(2)*w%r(2) + tz%r(2)*w%r(3)
            dir%r(3) =
     *          tx%r(3)*w%r(1) + ty%r(3)*w%r(2) + tz%r(3)*w%r(3)
         else
            dr%r(1) = 0.
            dr%r(2) = 0.
            dr%r(3) = 0.
            dir = atrack%vec%w
         endif
      else
         dr%r(1) = leng*atrack%vec%w%r(1)
         dr%r(2) = leng*atrack%vec%w%r(2)
         dr%r(3) = leng*atrack%vec%w%r(3)
         dir = atrack%vec%w
      endif

Here is the call graph for this function:

Here is the caller graph for this function: