C   04/03/80 704231727  MEMBER NAME  TFCTD1   (S)           FORTRAN
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
       SUBROUTINE TFCTD1(INPA)
C
C              COMPUTES THE HIT-POSITION AND PATH-LENGTH,AND
C          MAKES TABLES
C            TRK(5,N)    : Z, PHAI, P(MEV), PATH, ZS
C            ITRK(5,CNT#): CNT#, MULTIPLICITY, TRK#1, TRK#2, ...
C            ITRC(TRK#)  : CNT# <=> TRK#
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
CC
C
      COMMON/EVRUN/NRUN,IEV,ITYP,NEV,IMISS
      COMMON/TOFHIT/XTOF(50),YTOF(50),ZTOF(50)
       COMMON/CWORK/NR,RAW(5,42),NC,ICRT1(5,42),NTRK,ICRT2(50),TRK(5,50)
     - ,ITRC(50),NTC,ITRK(5,42),INFM(4),IR(14,50)
C   12/03/79 007281422  MEMBER NAME  DDATA    (S)           FORTRAN
C
      COMMON /BCS/ IDATA(10000)
      DIMENSION DATA(10000)
      INTEGER *2 HDATA(20000)
      EQUIVALENCE (IDATA(1),DATA(1)), (HDATA(1),IDATA(1))
C
       COMMON/CTPCNS/BENGEV,BKGAUS
C
       DATA PAI/ 3.1415927/ ,IBUG/100/,JSQRT/10/
       DATA RTOF,ZLN,DLMAX/ 920.0, 1610.0, 1.0E4/
       DATA DPAI/.14960/
      DATA IENTRY/0/
      IF(IENTRY.NE.0)  GOTO  49
      RTOF2= RTOF*RTOF
      IENTRY = 1
C
   49 IBUG = IBUG + 1
      IF(IBUG.LE.10) PRINT 105,NEV
  105 FORMAT(//' NEW EVENT TFCTD1 ',I5/)
      CALL SETSL(ITRC,0,50*4,0)
      CALL SETSL(TRK,0,250*4,0)
      CALL SETSL(ITRK,0,210*4,0)
C
      NTRK = IDATA(INPA+2)
      LTR = IDATA(INPA+3)
      ISTART = INPA+IDATA(INPA+1)
      IF(IBUG.LE.0) PRINT 676,NTRK,INPA,ISTART
  676 FORMAT(' TRACKS TFCTD1 ',3I10)
      IF(NTRK.GT.50) NTRK = 50
      IF(NTRK.LE.0) RETURN
C/////////
C
       DO 1000 I= 1,NTRK
C
C=========   LOAD TRACK INFORMATION
C
      IP=ISTART+(I-1)*LTR
            ITNUM= IDATA(IP+1)
C/////////
            IF(ITNUM.NE.I) GO TO 1000
C CHECK THAT TRACK HITS TOF COUNTER IN Z
      DZDR = DATA(IP+30)
      ZT   = DATA(IP+31) + DZDR * RTOF
      ZTOF(I) = ZT
      IF(ABS(ZT).GT.ZLN+100.)  GOTO  1000
            X1= DATA(IP+5)
            Y1= DATA(IP+6)
            Z1= DATA(IP+7)
            IF(ABS(Z1).GT.ZLN) GO TO 1000
            X2= DATA(IP+12)
            Y2= DATA(IP+13)
            Z2= DATA(IP+14)
            IF(ABS(Z2).GT.ZLN) GO TO 1000
C
C  COMPUTE CENTER OF CURVATURE FROM TRACK PARAMETERS
C
      IF(ABS(DATA(IP+25)).LT.1.E-28) GOTO  1000
      R = 1./DATA(IP+25)
      ITYPE = IDATA(IP+18)
      IF(IBUG.LE.10) PRINT 111,ITYPE,R
  101 FORMAT(1X,2I5,10F12.5)
  111 FORMAT(1X,I5,10F12.5)
      IF(IDATA(IP+18).EQ.2)  GOTO  30
      IF(IDATA(IP+18).NE.1)  GOTO  1000
C
C  CIRCLE
C
      IF(ABS(DATA(IP+19)).LT.1.E-28) GOTO  1000
      RC= 1./DATA(IP+19)
      D0 = DATA(IP+20)+RC
      ALPHC = DATA(IP+21)
      XM = D0*COS(ALPHC)
      YM = D0*SIN(ALPHC)
      XCL = DATA(IP+20)*COS(ALPHC)
      YCL = DATA(IP+20)*SIN(ALPHC)
      IF(IBUG.LE.0) PRINT 101,IDATA(IP+18),RC,D0,ALPHC,XM,YM,XCL,YCL
      GOTO  31
C
C PARABOLA
   30 CONTINUE
      RP = R
      ALPHP = DATA(IP+19)
      X0 = DATA(IP+20)
      Y0 = DATA(IP+21)
      XM = X0 + RP*SIN(ALPHP)
      YM = Y0 - RP*COS(ALPHP)
      D0 = SQRT(XM**2+YM**2)
      DS = D0-ABS(RP)
      IF(ABS(XM).LT..0000001)  GOTO  32
      ALPH = ATAN2(YM,XM)
      GOTO  33
   32 ALPH = PAI/2.
   33 XCL = DS*COS(ALPH)
      YCL = DS*SIN(ALPH)
      IF(IBUG.LE.0) PRINT 101,IDATA(IP+18),RP,ALPHP,X0,Y0,XM,YM,XCL,
     1 YCL
   31 CONTINUE
      ZCL = SQRT(XCL**2+YCL**2)*DZDR+DATA(IP+31)
C
C  COMPUTE HIT POINT AT TOF COUNTER =
C  POINT ON THE CIRCLE WITH RADIUS R, MIDDLE XM,YM WHICH HAS A DISTANC
C  FROM THE ORIGIN OF RTOF
C
C   STRAIGHT LINE EXTRAPOLATION
C
      RL = SQRT(X2*X2+Y2*Y2)
      RD = RTOF - RL
      IF(ABS(DATA(IP+15)).LT..00001) GOTO  650
      TANXY = DATA(IP+16)/DATA(IP+15)
      PHIXY = ABS(ATAN(TANXY))
      XTE = X2+RD*COS(PHIXY)*SIGN(1.,DATA(IP+15))
      YTE = Y2+RD*SIN(PHIXY)*SIGN(1.,DATA(IP+16))
      GOTO  651
  650 XTE = 0.
      YTE = RTOF
  651 IF(IBUG.LE.0) PRINT 102 ,XTE,YTE
C  FOR HIGH MOMENTA TAKE STRAIGHT EXTRAPOLATION
      IF(ABS(R).GT.7000.)  GOTO  50
C  FOR LOW MOMENTA CALCULATE INTERCEPT OF 2 CIRCLES
      A = RTOF2-R*R+XM*XM+YM*YM
      P = -YM*A/(XM**2+YM**2)
      Q = -(XM**2*RTOF2-A*A/4.)/(XM**2+YM**2)
      IF(IBUG.LE.0) PRINT 102,A,P,Q
  102 FORMAT(1X,10E12.5)
      IF(P*P-4.*Q.LE.0.)  GOTO  987
      YT1 = .5*(-P+SQRT(P*P-4.*Q))
      YT2 = .5*(-P-SQRT(P*P-4.*Q))
      DY1 = ABS(YTE-YT1)
      DY2 = ABS(YTE-YT2)
C TAKE SOLUTION WHICH IS CLOSER TO STRAIGHT EXTRAPOLATION
      DC11 = 5000.
      DC12 = 5000.
      IF(RTOF2-YT1*YT1.LT.0.) GOTO  7
      XT11= SQRT(RTOF2-YT1*YT1)
      XT12= -XT11
      DC11= (XT11-XM)**2+(YT1-YM)**2-R*R
      DC12= (XT12-XM)**2+(YT1-YM)**2-R*R
    7 IF(RTOF2-YT2*YT2.GT.0.) GOTO  8
      DC21 = 5000.
      DC22 = 5000.
      IF(DC11.GE.5000.)  GOTO  988
      GOTO  9
    8 XT21= SQRT(RTOF2-YT2*YT2)
      XT22= -XT21
      DC21= (XT21-XM)**2+(YT2-YM)**2-R*R
      DC22= (XT22-XM)**2+(YT2-YM)**2-R*R
    9 IF(IBUG.LE.0) PRINT 102,YT1,XT11,XT12,DC11,DC12,YT2,XT21,XT22,
     1 DC21,DC22
C DECIDE WHICH 2 SOLUTIONS FIT CIRCLES
      XT1 = XT11
      IF(ABS(DC12).LT.ABS(DC11)) XT1 = XT12
      XT2 = XT21
      IF(ABS(DC22).LT.ABS(DC21)) XT2 = XT22
C
C  TAKE SOLUTION CLOSEST TO STRAIGHT TRACK EXTRAPOLATION
C
      DE1 = SQRT((XTE-XT1)**2+(YTE-YT1)**2)
      DE2 = SQRT((XTE-XT2)**2+(YTE-YT2)**2)
      IF(IBUG.LE.0) PRINT 102,DE1,DE2
      YT = YT1
      IF(ABS(DE2).LT.ABS(DE1)) YT = YT2
      XT = XT1
      IF(ABS(DE2).LT.ABS(DE1)) XT = XT2
      XTOF(I) = XT
      YTOF(I) = YT
      IF(IBUG.LE.10) PRINT 102,XT,XTE,YT,YTE
      EXTR = SQRT((XT-XTE)**2+(YT-YTE)**2)
      GOTO 52
  987 CONTINUE
      ISQRT = ISQRT + 1
      IF(ISQRT.LE.0) PRINT 103,I,ISQRT
  103 FORMAT(' NEG SQRT ',5I10)
      IF(ISQRT.LE.0) PRINT 102,A,P,Q,R,X1,Y1,X2,Y2,XM,YM
      GOTO  1000
  988 JSQRT = JSQRT + 1
      IF(JSQRT.LE.0) PRINT 507,RTOF,YT
  507 FORMAT(' RTOF YT  --  NEG SQRT ',2E12.3)
      IF(JSQRT.LE.0) PRINT 102,A,P,Q
      IF(JSQRT.LE.0) PRINT 102,YT1,YT2,DY1,DY2,YT
      GOTO  1000
C
C=========   PATH LENGTH
C
   50 CONTINUE
      XT = XTE
      YT = YTE
      XTOF(I) = XT
      YTOF(I) = YT
   52 SEHNE = SQRT((XT-XCL)**2+(YT-YCL)**2)
      SINM = ABS(SEHNE/2./R)
      IF(SINM.GT.1.)  GOTO  1000
      ANGM = ASIN(SINM)
      PATHR = ABS(2.*ANGM*R)
      COSA = 1./SQRT(1.+DZDR*DZDR)
      IF(COSA.EQ.0.) COSA = 1.
      PATH= PATHR/COSA
      PHAI = ATAN2(YT,XT)
      IF(PHAI.LT.0.) PHAI = 2.*PAI+PHAI
      IF(IBUG.LE.0) PRINT 102,SEHNE,PATHR,PATH,PHAI
C/////////
C
C=========   STORE THE RESULTS INTO COMMON
C---   ITRC(50) : TRACK # I <=> COUNTER # NCNT
            PHD= PHAI-DPAI/2.
            IF(PHD.GT.0.) GO TO 300
            NCNT= 1
            GO TO 310
 300        NCNT= IFIX(PHD/DPAI)+2
            IF(NCNT.GE.43) NCNT=NCNT-42
 310        ITRC(I)= NCNT
      IF(IBUG.LE.10) PRINT 101,I,NCNT,ZCL
      IF(IBUG.LE.10) PRINT 111,I,XTOF(I),YTOF(I),ZTOF(I)
C---   ITRK(5,42) : COUNTER # NCNT <=> TRACK MULTIPLICITY
            ITRK(1,NCNT)= NCNT
            MTPY= ITRK(2,NCNT)
            MTPY= MTPY+1
            ITRK(2,NCNT)= MTPY
       IF(MTPY.GE.4) GO TO 320
            MTPY= MTPY+2
            ITRK(MTPY,NCNT)= I
C---   TRK(5,50) : TRACK # <=> Z,PHAI,PATH,MOMENTUM,Z-SCATT
C THE MINUS SIGN IS TO COMPENSATE FOR THE DEFINITION OF R ABOVE
 320        PGEVR= -0.03*R*ABS(BKGAUS)/1000.
      PGEV = PGEVR/COSA
      IF(IBUG.LE.10) PRINT 751,R,PGEVR,COSA,PGEV
  751 FORMAT( ' PGEV ',4E12.5/ )
            TRK(1,I)= ZTOF(I)
            TRK(2,I)= PHAI
            TRK(3,I)= PATH
            TRK(4,I)= PGEV
            TRK(5,I)= ZCLOS
C
C/////////
 1000   CONTINUE
C/////////
        RETURN
        END
