C These routines calculate resolutions in various track parameters due
C to both position resolution and multiple scattering individually as
C well as the total (i. e., their sum in quadrature).
C
C The subroutines are:
C
C      SUBROUTINE REZEST_CURVE(DET, LAMBDA, L, EPSILON, N_M, P, M, N_RL,
C     X     B, DK_MS, DK_RES, DK)
C
C      SUBROUTINE REZEST_AZIMUTH(DET, LAMBDA, L, EPSILON, N_M, P, M,
C     X     N_RL, N_RL_FRONT, B, R_M, DK, DPHI_MS, DPHI_RES, DPHI_K,
C     x     DPHI)
C
C      SUBROUTINE REZEST_POLAR(DET, LAMBDA, L, EPSILON, STEREO,
C     X     N_M, P, M, N_RL, N_RL_FRONT, DTHETA_MS, DTHETA_RES,
C     X     DTHETA)
C
C All have the following REAL*4 input arguments in common:
C
C     DET      Detector, CHARACTER*3, either 'CDC' or 'FDC'
C     LAMBDA   Dip angle, difference in polar angle in lab between track
C              and pi/2 (i. e., 90 degrees) (radians)
C     L_DET    Length of detector
C              CDC: the distance from the inner radius of the chamber to the
C                   outer radius,
C              FDC: the distance from the front of the chamber to the back of
C                   the chamber
C              (meters)
C     EPSILON  Position resolution (meters)
C     N_M      Number of measurements (real number)
C     P        Magnitude of total momentum (GeV/c)
C     M        Mass of the particle (GeV)
C     N_RL     Number of radiation lengths in the detector (real number)
C
C In addition REZEST_CURVE and REZEST_AZIMUTH require as input the REAL*4
C argument:
C
C     B        Magnetic field (Tesla)
C
C REZEST_AZIMUTH and REZEST_POLAR require a REAL*4 argument:
C
C     N_RL_FRONT Number of radiations lengths in front of the detector
C
C Finally, REZEST_AZIMUTH requires a REAL*4 argument:
C
C     DK       Resolution in curvature, as output from REZEST_CURVE (see below)
C              (inverse meters)
C
C Output arguments are all REAL*4:
C
C REZEST_CURVE: DK_MS, DK_RES, DK
C               curvature resolution due to multiple scattering, position
C               resolution and their sum in quadrature respectively
C               (inverse meters)
C
C REZEST_AZIMUTH: DPHI_MS, DPHI_RES, DPHI
C               resolution in azimuthal angle due to multiple scattering,
C               position resolution and their sum in quadrature respectively
C               (radians)
C
C REZEST_POLAR: DTHETA_MS, DTHETA_RES, DTHETA
C               resolution in polar angle due to multiple scattering,
C               position resolution and their sum in quadrature respectively
C               (radians)
C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
      SUBROUTINE REZEST_CURVE(DET, LAMBDA, L_DET, EPSILON, N_M, P, M,
     X     N_RL, B, DK_MS, DK_RES, DK)
C     Calculates resolution in curvature
      IMPLICIT NONE
      INCLUDE "rezest.par"
      CHARACTER*3 DET
      REAL*4 LAMBDA, L_DET, EPSILON, N_M, P, M, N_RL, B
      REAL*4 E, GAMMA, BETA
      REAL*4 L, L_PRIME, DK_RES, DK_MS, DK, N_RL_TOTAL

      IF (N_M .LE. N_M_MIN) THEN
         DK_MS = BIG_VALUE
         DK_RES = BIG_VALUE
         DK = BIG_VALUE
         RETURN
      END IF

      CALL KINEMATICS(P, M, E, GAMMA, BETA)
      CALL SET_LENGTHS(DET, L_DET, LAMBDA, L, L_PRIME)
      dk_res = epsilon/L_prime**2*sqrt(720/(n_m + 4))
      IF (DET .EQ. 'CDC') THEN
         N_RL_TOTAL = N_RL/COS(LAMBDA)
      ELSE IF (DET .EQ. 'FDC') THEN
         N_RL_TOTAL = N_RL/COS(3.14159/2.0 - LAMBDA)
      ELSE
         WRITE (6,*) 'BAD DETECTOR TYPE'
         STOP   
      ENDIF
      dk_ms = 0.016/(L*p*beta*cos(lambda)**2)*sqrt(n_rl_TOTAL)
      dk = sqrt(dk_res**2 + dk_ms**2)
      RETURN
      END
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
      SUBROUTINE REZEST_AZIMUTH(DET, LAMBDA, L_DET, EPSILON, N_M, P, M,
     X     N_RL, N_RL_FRONT, B, R_M, DK, DPHI_MS, DPHI_RES, DPHI_K,
     x     DPHI)
      IMPLICIT NONE
      INCLUDE "rezest.par"
      CHARACTER*3 DET
      REAL*4 LAMBDA, L_DET, EPSILON, N_M, P, M, N_RL, N_RL_FRONT, B,
     X     R_M, DK
      REAL*4 E, GAMMA, BETA
      REAL*4 THETA0, C, DPHI_MS, DPHI_RES, DPHI_K, DPHI,
     X     N_RL_TOTAL, L, L_PRIME, PT, K, ALPHA
      DATA C/1.0/

      IF (N_M .LE. N_M_MIN) THEN
         DPHI_MS = BIG_VALUE
         DPHI_RES = BIG_VALUE
         DPHI_K = BIG_VALUE
         DPHI = BIG_VALUE
         RETURN
      END IF

      CALL SET_LENGTHS(DET, L_DET, LAMBDA, L, L_PRIME)
      CALL KINEMATICS(P, M, E, GAMMA, BETA)
      IF (DET .EQ. 'CDC') THEN
         N_RL_TOTAL = (N_RL + N_RL_FRONT)/COS(LAMBDA)
      ELSE IF (DET .EQ. 'FDC') THEN
         N_RL_TOTAL = (N_RL + N_RL_FRONT)/SIN(LAMBDA)
      ELSE
         WRITE (6,*) 'BAD DETECTOR TYPE'
         STOP   
      ENDIF
      THETA0 = 0.0136/(BETA*C*P)*SQRT(N_RL_TOTAL)
     X     *(1.0 + 0.038*LOG(N_RL_TOTAL))
      DPHI_MS = THETA0/SQRT(3.0)
      IF (N_M .GE. 2.0) THEN
         DPHI_RES = EPSILON/L_PRIME*SQRT(12*(N_M - 1)/(N_M*(N_M + 1)))
      ELSE
         DPHI_RES = BIG_VALUE
      END IF
      PT = P*COS(LAMBDA)
      K = 0.3*B/PT
      IF (2.0/K .LE. R_M) THEN
C        linear track approximation breakdown: 2 X radius of curvature less
C        than measurement radius
         DPHI_K = BIG_VALUE
      ELSE
         ALPHA = ASIN(0.5*R_M*K)
         DPHI_K = R_M*DK/COS(0.5*ALPHA)
      ENDIF
      DPHI = SQRT(DPHI_RES**2 + DPHI_MS**2 + DPHI_K**2)
      RETURN
      END
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
      SUBROUTINE REZEST_POLAR(DET, LAMBDA, L_DET, EPSILON, STEREO,
     X     N_M, P, M, N_RL, N_RL_FRONT, DTHETA_MS, DTHETA_RES,
     X     DTHETA)
      IMPLICIT NONE
      INCLUDE "rezest.par"
      CHARACTER*3 DET
      REAL*4 LAMBDA, L_DET, EPSILON, STEREO, N_M, P, M, N_RL,
     X     N_RL_FRONT
      REAL*4 E, GAMMA, BETA
      REAL*4 C, N_RL_TOTAL, THETA0, EPS_EFF, DTHETA_MS, DTHETA_RES,
     x     DTHETA, PIOVER2

      DATA C/1.0/
      DATA PIOVER2 /1.5707963/

      IF (N_M .LE. N_M_MIN) THEN
         DTHETA_MS = BIG_VALUE
         DTHETA_RES = BIG_VALUE
         DTHETA = BIG_VALUE
         RETURN
      END IF

      CALL KINEMATICS(P, M, E, GAMMA, BETA)
      IF (DET .EQ. 'CDC') THEN
         N_RL_TOTAL = (N_RL + N_RL_FRONT)/COS(LAMBDA)
      ELSE IF (DET .EQ. 'FDC') THEN
         N_RL_TOTAL = (N_RL + N_RL_FRONT)/COS(PIOVER2 - LAMBDA)
      ELSE
         WRITE (6,*) 'BAD DETECTOR TYPE'
         STOP   
      ENDIF
      THETA0 = 0.0136/(BETA*C*P)*SQRT(N_RL_TOTAL)
     X     *(1.0 + 0.038*LOG(N_RL_TOTAL))
      DTHETA_MS = THETA0/SQRT(3.0)
      IF (DET .EQ. 'CDC') THEN
         EPS_EFF = EPSILON*COS(LAMBDA)**2/SIN(STEREO)
      ELSE IF (DET .EQ. 'FDC') THEN
         EPS_EFF = EPSILON*COS(PIOVER2 - LAMBDA)**2
      ELSE
         WRITE (6,*) 'BAD DETECTOR TYPE'
         STOP
      ENDIF
      IF (N_M .GE. 2.0) THEN
         DTHETA_RES = EPS_EFF/L_DET
     X        *SQRT(12*(N_M - 1)/(N_M*(N_M + 1)))
      ELSE
         DTHETA_RES = BIG_VALUE
      END IF
      DTHETA = SQRT(DTHETA_RES**2 + DTHETA_MS**2)
      RETURN
      END
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
      SUBROUTINE KINEMATICS(P, M, E, GAMMA, BETA)
      REAL*4 P, M
      REAL*4 E, GAMMA, BETA
      REAL*4 C
      DATA C/1.0/
      E = sqrt(p**2 + m**2)
      gamma = E/m
      beta = sqrt(1.0 - 1/gamma**2)
      RETURN
      END
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
      SUBROUTINE SET_LENGTHS(DET, L_DET, LAMBDA, L, L_PRIME)
      REAL*4 L_DET, LAMBDA, L, L_PRIME
      CHARACTER*3 DET
      IF (DET .EQ. 'CDC') THEN
         L_PRIME = L_DET
         L = L_DET/COS(LAMBDA)
      ELSE IF (DET .EQ. 'FDC') THEN
         L_PRIME = L_DET/TAN(LAMBDA)
         L = L_DET/SIN(LAMBDA)
      ELSE
         STOP 'BAD DETECTOR'
      END IF
      RETURN
      END
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC