*--------1---------2---------3---------4---------5---------6---------7--

* This module tabulates the track curvature
* in the solenoidal magnetic field

      subroutine HELIX

      include './Include/padgeom.inc'
      include './Include/trkdata.inc'

      double precision PT, RanTheta, RanPhi
      double precision Vlong, Vtrans, Rbend, Aspeed
      double precision Xspacal, Yspacal, Tspacal

      double precision Todd(4), Teven(4)
      double precision EvenX(4), OddX(4), EvenY(4), OddY(4)

      double precision Slight, H1field, Pi

      parameter (Pi = 3.14159265358979)

* Speed of light mm/sec
*----------------------------------
      parameter (Slight = 299792458000.)

* Magnetic field strength, T
* (vector B is collinear with Z)
*----------------------------------
      parameter (H1field = 1.15)

* For relativistic electrons
* (E >> 511 keV) the statement E=pc
* is valid to a very good precision,
* therefore V=c
*----------------------------------
      PT = sqrt(PX**2 + PY**2)

* Z-axis is inverted for simplicity in
* calculations, thus a complementary to
* Theta angle is being considered
*----------------------------------
      RanTheta = Pi - atan2(PT, PZ)
      RanPhi = atan2(PY, PX)

      Vlong = Slight * cos(RanTheta)
      Vtrans = Slight * sin(RanTheta)

      Rbend = 1000. * PT / 0.3 / H1field
      Aspeed = Vtrans / Rbend

* Sign of the track curvature
* for electrons or positrons
*----------------------------------
      Aspeed = -1. * Aspeed * sign(1, Pcode)
      RanPhi = RanPhi + Pi / 2. * sign(1, Pcode)

* Track projection onto the SPACAL plane:
*----------------------------------
      Tspacal = (SpacalPos + ZV) / Vlong
      Xspacal = XV + Rbend*(cos(Aspeed*Tspacal + RanPhi) - cos(RanPhi))
      Yspacal = YV + Rbend*(sin(Aspeed*Tspacal + RanPhi) - sin(RanPhi))
      Rspacal = real(SQRT(Xspacal**2 + Yspacal**2))

* Tabulating X,Y coordinates in every Z-plane:
*----------------------------------
      do I = 1, 4
        Teven(I) = (EvenPlanePos(I) + ZV) / Vlong
        Todd(I) = (OddPlanePos(I) + ZV) / Vlong

        OddX(I)  = XV + Rbend*(cos(Aspeed*Todd(I) +RanPhi)-cos(RanPhi))
        EvenX(I) = XV + Rbend*(cos(Aspeed*Teven(I)+RanPhi)-cos(RanPhi))

        OddY(I)  = YV + Rbend*(sin(Aspeed*Todd(I) +RanPhi)-sin(RanPhi))
        EvenY(I) = YV + Rbend*(sin(Aspeed*Teven(I)+RanPhi)-sin(RanPhi))

        EvenR(I) = SQRT(EvenX(I)**2 + EvenY(I)**2)
        OddR(I) = SQRT(OddX(I)**2 + OddY(I)**2)

        EvenPhi(I) = 180. * atan2(EvenY(I), EvenX(I)) / Pi
        OddPhi(I) = 180. * atan2(OddY(I), OddX(I)) / Pi

        if (EvenPhi(I).lt.0) EvenPhi(I) = EvenPhi(I) + 360.
        if (OddPhi(I).lt.0) OddPhi(I) = OddPhi(I) + 360.

      end do

      return
      end