      SUBROUTINE COSMIC_UTC
C
C  This is a dplot to look at cosmic ray data.
C
C  15-Feb-1996   Modified the dphi vs tanl1 plot so that we can look for
C                 a twist in the anode wires.
C                Use PARAM(1) to disable RD_TRK (gives thousands of errors
C                 when TD's not read out).
C  24-Feb-1996   Modified the way we extrapolate the track to the foils
C                 when plotting z residual vs strip number.  Previously
C                 we used TRKPOS, which has a problem with 'direction'
C                 for cosmics fit to a single track.  We now use the
C                 turning angle.
C  18-Mar-1996   Largely restructured.  We now do both single-helix fits
C                 and independent track fits, depending on what quantity
C                 we want to look at.
C  22-Mar-1996   Added sector-by-sector RSSC z offsets
C  29-Mar-1996   Added extrapolation to B4
C  05-Apr-1996   Fixed nhitpl calculation for 'cosmic' option.
C  22-Apr-1997   Added stuff to check the foil stereo angles
C  05-Dec-1997   Added stuff to look at possible twist and tilt of SL 2

C-- Global declarations

      IMPLICIT NONE
#include <fail.par>
#include <ss.par>
#include <dc_track.par>
#include <luns.cmn>
#include <info.cmn>
#include <hcuts.cmn>
#include <dc_track.cmn>
#include <cfail.cmn>
#include <utc_geom.cmn>
#include <utc_clusters.cmn>    ! has to precede utc_candidate
#include <utc_candidate.cmn>
#include <utc_tdcdata.cmn>
#include <utc_adcdata.cmn>
#include <utc_xyz.cmn>
#include <utc_wires.cmn>
#include <utc_strips.cmn>
#include <utc_print.cmn>
#include <utc_xy_pat.cmn>
#include <utc_xy_cuts.cmn>
#include <utc_z_cuts.cmn>
#include <ss_index.cmn>
#include <fail.cmn>
#include <rdtrack.cmn>
#include <rd_anal.cmn>
#include <adcamp.cmn>
#include <tr2bits.cmn>

C-- HBOOK/PAW stuff 

      INTEGER*4 HBSIZE
      PARAMETER(HBSIZE=5000000)
      COMMON/PAWC/HMEMOR(HBSIZE)
      REAL*4 HMEMOR

      INTEGER*4 NPRIME
      PARAMETER(NPRIME=1024)

      INTEGER*4 NUMT
      PARAMETER(NUMT=2)
      REAL*4    XFIL(NUMT)
      CHARACTER*8 NAMES(NUMT)

      INTEGER*4 IOFF,IOFF2,IOFF3
      PARAMETER(IOFF=1000,IOFF2=2000,IOFF3=3000)

C-- Other local declarations

      INTEGER*4 MNHPT
      PARAMETER(MNHPT=24)                 ! Max number of hits per track
      INTEGER*4 MNHPL
      PARAMETER(MNHPL=2)                  ! Max number of hits per layer

      REAL*4    RIC
      PARAMETER(RIC=6.32)                 ! Radius of middle of I-counter
      REAL*4    RTC
      REAL*4    RBV
      REAL*4    RB4
      
      INTEGER*4 IERR
      INTEGER*4 IC,IP,IPTDC
      INTEGER*4 NH,NHIT,NHITPL(24)
      INTEGER*4 IHIT
      INTEGER*4 ICELL,ILAY,ISLAY,NLAYER,JLAYER
      INTEGER*4 IWIRE(MNHPT)
      INTEGER*4 I,J
      INTEGER*4 NDOF
      INTEGER*4 ITRACK,ITRK1,ITRK2
      INTEGER*4 IS,IND,ISTRIP,IFOIL,IW
      INTEGER*4 IHITPL(MNHPL,24)
      INTEGER*4 NLAYFIT,ILAYER
      INTEGER*4 ISTRIP_MIN,ICELL_MIN
      INTEGER*4 ISEC1,ISEC2,ISECT(4),ISECIC,ISECTC,ISECBV,ISEC,IEND
      INTEGER*4 IDUMMY
      INTEGER*4 NHITZ_MIN
      INTEGER*4 IMOD,JSEC
      INTEGER*4 IZBIN,ILR

      REAL*4    R,CHIPXY,CHIDFXY,PHI0,CHIXY,CHIDFZ,CHIPZ
      REAL*4    TWOPI,PI,PIBY2
      REAL*4    TIME(MNHPT)
      REAL*4    TRKPAR(3),XYPCA(2),DOCA,PHIVU,PHIDCA
      REAL*4    RESID,DDRIFT,RESID2
      REAL*4    PNEW(3),ERRPAR(6),CHIDFNEW
      REAL*4    TANL,TANLT,SINL,COSL
      REAL*4    MYRHO,MYD0,MYPHI0
      REAL*4    XT,YT,RT,PHIT,DHT
      REAL*4    DR,DRMIN,DP,DPMIN
      REAL*4    XTOLD,YTOLD
      REAL*4    XYPCA2(2,MNHPT),DOCA2(MNHPT)
      REAL*4    PHIVU2(MNHPT),PHIDCA2(MNHPT)
      REAL*4    PHW,PHS
      REAL*4    XYRESMAX,ZRESMAX
      REAL*4    XFOIL,YFOIL,ZFOIL,DXFOIL,DYFOIL,DZFOIL
      REAL*4    TMIN,DZMIN,DZ
      REAL*4    SINL1,COSL1,SINL2,COSL2,TANL1,TANL2
      REAL*4    DELTAP,DELTAD
      REAL*4    XIC,YIC,ZIC,COSIC(3),PHIIC,EIC
      REAL*4    XTC,YTC,ZTC,COSTC(3),PHITC,ETC
      REAL*4    XB4,YB4,ZB4(2),COSB4(3),PHIB4,EB4
      REAL*4    XBV,YBV,ZBV,COSBV(3),PHIBV,EBV,EBVLAST
      REAL*4    PHI1,PHI2,D01,D02,Z01,Z02
      REAL*4    DUMMY
      REAL*4    RSC1,RSC2,XSC1,XSC2,YSC1,YSC2,ZSC1,ZSC2
      REAL*4    DXSC1,DYSC1,DZSC1,DXSC2,DYSC2,DZSC2
      REAL*4    PHIUT1,PHIUT2,PHISC1,PHISC2
      REAL*4    DPHI1,DPHI2,DUTSC1,DUTSC2,DZ1,DZ2
      REAL*4    Z1,Z2
      REAL*4    PATH,RSC(4),DPHI(4),ZED(4),PARAMS(2),COV(2,2),ZRESID(4)
      REAL*4    XTRAER(4)
      REAL*4    CHIZ,SINP,COSP,DENOM,SINTH,COSTH
      REAL*4    PHIW1,DPHIW
      REAL*4    TLMAX
      REAL*4    RSSC_ZOFF(24,2)
      REAL*4    R_LAY(12)
      REAL*4    DX,DY,X,Y,Z
    
      LOGICAL*4 LHIT(MNHPT)
      LOGICAL*4 LMC
      LOGICAL*4 LFIT(4)
      LOGICAL*4 LFUP
      LOGICAL*4 LFIRST
      REAL PROB_XY,PROB_Z
      CHARACTER*2 CNUM
      
C-- Parameters

      INTEGER*4 LDST
      PARAMETER(LDST=99)

C-- Functions

      REAL*4    PROB

C-- Data statements

      DATA LFIRST/.TRUE./
      DATA LFIT/4*.TRUE./
      DATA RSSC_ZOFF/-1.3, -1.4, -0.8, -1.3, -1.1, -0.8,
     &               -0.8, -0.7, -1.4, -0.9, -1.4, -1.6,
     &               -1.0, -1.8, -1.1, -1.4, -1.2, -1.1,
     &               -1.5, -1.0, -0.9, -1.0, -1.2, -1.2,
     &               -1.0, -1.5, -1.3, -0.9, -0.8, -1.0,
     &               -1.0, -1.0, -0.8, -1.1, -0.9, -1.3,
     &               -1.4, -1.1, -0.8, -1.0, -1.5, -0.8,
     &               -0.4, -1.0, -1.8, -1.2, -1.7, -0.9/

      SAVE

C ***************************************
C  Global init entrypoint
C ***************************************

      ENTRY DEFINE

C-- Initialize some constants

      PIBY2 = ASIN(1.0)
      PI    = 2.0*PIBY2
      TWOPI = 2.0*PI

      IF(UTCXYRES.NE.0.) THEN
        XYRESMAX = 10.*UTCXYRES
      ELSE
        XYRESMAX = 0.1
      ENDIF

      IF(UTCZRES.NE.0.) THEN
        ZRESMAX = 10.*UTCZRES
      ELSE
        ZRESMAX = 3.0
      ENDIF

C-- Number of foil hits required for z residual plots

      NHITZ_MIN = PARAM(2)
      
C-- Initialize HBOOK

      CALL HLIMIT(HBSIZE)

C-- Open the HBOOK RZ file given by KOFIA_HIS

      CALL HBKINI(IERR)
      IF(IERR.NE.0) THEN
        CALL KERROR(4,1,'define','Problem opening hbook file!!!')
        STOP
      ENDIF

C-- Book the ntuples

c-- Book the histograms

C-- Single-helix fit
      CALL HBOOK1(IOFF+1,'Number of tracks per event',10,-0.5,9.5,0.)
      CALL HBOOK1(IOFF+2,'Number of xy hits per track',30,-0.5,29.5,0.)
      CALL HBOOK1(IOFF+3,'Number of z hits per track',20,-0.5,19.5,0.)
      CALL HBOOK1(IOFF+4,'Number of xy layers per track',
     &                                             30,-0.5,29.5,0.)
      CALL HBOOK1(IOFF+5,'Track charge',4,-2.5,1.5,0.)
      CALL HBOOK1(IOFF+6,'XY chisq prob',100,0.,1.,0.)
      CALL HBOOK1(IOFF+7,'XY Chisq/ndof',100,0.,20.,0.)
      CALL HBOOK1(IOFF+8,'Z chisq prob',100,0.,1.,0.)
      CALL HBOOK1(IOFF+9,'Z Chisq/ndof',100,0.,20.,0.)
      CALL HBOOK2(IOFF+10,'Impact parameter vs phi0',100,0.,6.28,
     &                    100,-30.,30.,0.)
      CALL HBOOK1(IOFF+11,'Track Z0',100,-50.,50.,0.)
      CALL HBOOK2(IOFF+12,'Impact parameter vs Z0',
     &             20,-20.,20.,50,-30.,30.,0.)

      CALL HBOOK1(IOFF+21,'Z chisq/ndof. Nz .ge. nhitzmin',
     &                                         100,0.,20.,0.)
      CALL HBOOK1(IOFF+22,'IC energy',100,0.,5.,0.)
      CALL HBOOK1(IOFF+23,'Z at I-counter',200,-25.,25.,0.)
      CALL HBOOK1(IOFF+24,'TC energy',100,0.,10.,0.)
      CALL HBOOK1(IOFF+25,'Z at T-counter',200,-50.,50.,0.)      
      CALL HBOOK1(IOFF+26,'B4 energy',100,0.,2.,0.)
      CALL HBOOK1(IOFF+27,'Z at B4. B4 hits',200,-25.,25.,0.)
      CALL HBOOK1(IOFF+28,'Z at B4. B4 misses',200,-25.,25.,0.)
      CALL HBOOK2(IOFF+101,'XY residuals (cm) by layer',
     &         100,-XYRESMAX,XYRESMAX,12,0.5,12.5,0.)
      CALL HBOOK2(IOFF+102,'Scaled XY residuals (cm) by layer',
     &         100,-10.,10.,12,0.5,12.5,0.)
      CALL HBOOK2(IOFF+103,'Pulse area vs XY residual (cm)',
     &         100,-XYRESMAX,XYRESMAX,100,0.5,5000.5,0.)
      CALL HBOOK2(IOFF+104,'Wire layer vs pulse area',
     &         100,0.5,5000.5,12,0.5,12.5,0.)
      CALL HBOOK2(IOFF+105,'No. xy hits vs XY residual (cm)',
     &         100,-XYRESMAX,XYRESMAX,25,5.5,30.5,0.)
      CALL HBOOK2(IOFF+111,'Z residuals (cm) by foil',
     &         100,-ZRESMAX,ZRESMAX,6,0.5,6.5,0.)
      CALL HBOOK2(IOFF+112,'Pulse area vs Z residual (cm)',
     &         100,-ZRESMAX,ZRESMAX,100,0.5,10000.5,0.)
      CALL HBOOK2(IOFF+113,'Foil vs pulse area',
     &         100,0.5,10000.5,6,0.5,6.5,0.)
      CALL HBOOK2(IOFF+114,'No. foils hit vs Z residual (cm)',
     &         100,-ZRESMAX,ZRESMAX,12,0.5,12.5,0.)

      DO IFOIL=1,6
        WRITE(CNUM,100) IFOIL
100     FORMAT(I2)
        CALL HBOOK2(IOFF+120+IFOIL,
     &      'Strip cluster vs wire pulse area. Foil'//CNUM,
     &         50,0.5,10000.5,50,0.5,10000.5,0.)
      END DO

      DO IFOIL=1,6
        IF(MOD(IFOIL,2).EQ.0) THEN
          ISLAY = IFOIL-1
        ELSE
          ISLAY = IFOIL
        ENDIF
        WRITE(CNUM,100) IFOIL
        CALL HBOOK2(IOFF+130+IFOIL,
     &      'Z vs Z residual (cm). Nhitz.ge.nhitzmin. Foil'//CNUM,
     &      100,-ZRESMAX,ZRESMAX,
     &      50,-0.5*ZLENL(ISLAY)-5.,0.5*ZLENL(ISLAY)+5.,0.)
        CALL HBOOK2(IOFF+140+IFOIL,
     &      'Strip vs Z residual (cm). Nhitz.ge.nhitzmin. Foil '//CNUM,
     &      100,-ZRESMAX,ZRESMAX,
     &      NSFOIL(IFOIL),0.5,FLOAT(NSFOIL(IFOIL))+0.5,0.)
      END DO
      
      CALL HBOOK2(IOFF+201,'Signed DOCA residuals (cm) by layer',
     &         100,-0.1,0.1,12,0.5,12.5,0.)
      CALL HBOOK1(IOFF+202,'DOCA residual - std. residual',
     &         100,-0.1,0.1,0.)

C-- Two branches fit separately
      
      CALL HBOOK1(IOFF2+1,'Number of tracks per event',10,-0.5,9.5,0.)
      CALL HBOOK1(IOFF2+2,'Number of RS tracks',10,-0.5,9.5,0.)
      CALL HBOOK1(IOFF2+3,'Number of xy hits per track',30,-0.5,29.5,0.)
      CALL HBOOK1(IOFF2+4,'Number of z hits per track',20,-0.5,19.5,0.)
      
C-- UTC-RSSC matching

      CALL HBOOK2(IOFF2+301,'Delta(xy). UTC-RSSC inner. By sector',
     &            100,-5.,5.,24,0.5,24.5,0.)
      CALL HBOOK2(IOFF2+302,'Delta(xy). UTC-RSSC outer. By sector',
     &            100,-5.,5.,24,0.5,24.5,0.)
      CALL HBOOK2(IOFF2+303,'Delta(z). UTC-RSSC inner. By sector',
     &            100,-20.,20.,24,0.5,24.5,0.)
      CALL HBOOK2(IOFF2+304,'Delta(z). UTC-RSSC outer. By sector',
     &            100,-20.,20.,24,0.5,24.5,0.)
      CALL HBOOK2(IOFF2+305,'Zutc vs Delta(z) UTC-RSSC inner',
     &            100,-20.,20.,120,-60.,60.,0.)
      CALL HBOOK2(IOFF2+306,'Zutc vs Delta(z) UTC-RSSC outer',
     &            100,-20.,20.,120,-60.,60.,0.)

C-- Correlations between the two branches
      
      CALL HBOOK2(IOFF2+211,'Tanl1 vs (tanl1-tanl2)',
     &            100,-0.5,0.5,50,0.,1.,0.)
      CALL HBOOK2(IOFF2+212,'Z01 vs (Z01-Z02)',
     &            100,-5.,5.,50,-25.,25.,0.)
      CALL HBOOK2(IOFF2+213,'Phi1 vs (phi1-phi2)',
     &            100,-0.05,0.05,50,0.,6.28,0.)
      CALL HBOOK2(IOFF2+214,'D01 vs (D01+D02)',
     &            100,-0.3,0.3,50,-10.,10.,0.)
      CALL HBOOK2(IOFF2+215,'Tanl1 vs (Z01-Z02)',
     &            100,-5.,5.,50,0.,1.,0.)
      CALL HBOOK2(IOFF2+216,'DZFoil6 vs (phi1-phi2)',
     &            100,-0.05,0.05,50,0.,50.,0.)
      CALL HBOOK2(IOFF2+217,'DZFoil6 vs (D01+D02)',
     &            100,-0.2,0.2,50,0.,50.,0.)
      CALL HBOOK2(IOFF2+218,
     &            'Tanl1 vs (tanl1-tanl2). Same-side cosmics',
     &            100,-0.5,0.5,50,0.,1.,0.)
      CALL HBOOK2(IOFF2+219,
     &            'Tanl1 vs (tanl1-tanl2). Opposite-side cosmics',
     &            100,-0.5,0.5,50,0.,1.,0.)

C-- Residuals in SL2 to look for tilt, twist

      CALL HBOOK2(IOFF2+101,'Signed DOCA residual vs Z. SL 2',
     &            100,-25.,25.,100,-0.2,0.2,0.)
      CALL HBOOK2(IOFF2+102,'Signed DOCA residual vs Z. Layer 5',
     &            100,-25.,25.,100,-0.2,0.2,0.)
      CALL HBOOK2(IOFF2+103,'Signed DOCA residual vs Z. Layer 6',
     &            100,-25.,25.,100,-0.2,0.2,0.)
      CALL HBOOK2(IOFF2+104,'Signed DOCA residual vs Z. Layer 7',
     &            100,-25.,25.,100,-0.2,0.2,0.)
      CALL HBOOK2(IOFF2+105,'Signed DOCA residual vs Z. Layer 8',
     &            100,-25.,25.,100,-0.2,0.2,0.)

      CALL HBOOK2(IOFF2+110,
     & 'Signed DOCA residual vs phi0. SL 2. Z=-25 to -20',
     & 100,0.,6.28,100,-0.2,0.2,0.)
      CALL HBOOK2(IOFF2+111,
     & 'Signed DOCA residual vs phi0. SL 2. Z=-20 to -15',
     & 100,0.,6.28,100,-0.2,0.2,0.)
      CALL HBOOK2(IOFF2+112,
     & 'Signed DOCA residual vs phi0. SL 2. Z=-15 to -10',
     & 100,0.,6.28,100,-0.2,0.2,0.)
      CALL HBOOK2(IOFF2+113,
     & 'Signed DOCA residual vs phi0. SL 2. Z=-10 to -5',
     & 100,0.,6.28,100,-0.2,0.2,0.)
      CALL HBOOK2(IOFF2+114,
     & 'Signed DOCA residual vs phi0. SL 2. Z=-5 to 0',
     & 100,0.,6.28,100,-0.2,0.2,0.)
      CALL HBOOK2(IOFF2+115,
     & 'Signed DOCA residual vs phi0. SL 2. Z=0 to 5',
     & 100,0.,6.28,100,-0.2,0.2,0.)
      CALL HBOOK2(IOFF2+116,
     & 'Signed DOCA residual vs phi0. SL 2. Z=5 to 10',
     & 100,0.,6.28,100,-0.2,0.2,0.)
      CALL HBOOK2(IOFF2+117,
     & 'Signed DOCA residual vs phi0. SL 2. Z=10 to 15',
     & 100,0.,6.28,100,-0.2,0.2,0.)
      CALL HBOOK2(IOFF2+118,
     & 'Signed DOCA residual vs phi0. SL 2. Z=15 to 20',
     & 100,0.,6.28,100,-0.2,0.2,0.)
      CALL HBOOK2(IOFF2+119,
     & 'Signed DOCA residual vs phi0. SL 2. Z=20 to 25',
     & 100,0.,6.28,100,-0.2,0.2,0.)
      
      RETURN

C ***************************************
C  Event analysis
C ***************************************

      ENTRY DPLOT

      IF(LFIRST) THEN
        LFIRST = .FALSE.
        DO ISLAY=1,3
          DO ILAY=1,4
            NLAYER = (ISLAY-1)*4 + ILAY
            R_LAY(NLAYER) = RLAY(ILAY,ISLAY)
          END DO
        END DO
      ENDIF
      
      CALL BATCH_LOG
      ICFAIL=0

      CALL TR2UNP(0)
      IF(BD_ECXE) THEN
        ICFAIL=1
        RETURN
      ENDIF
      IF(.NOT.BD_COSM)RETURN
C-- First fit the two branches of the cosmic ray to a single track

      COSMIC=.TRUE.
      CALL UTC_TRACK(.TRUE.,.TRUE.)
      
      CALL HFILL(IOFF+1,FLOAT(NTRACK_D),0.,1.)

      IF(IFAIL_CODE(I_UT).NE.0) THEN
        ICFAIL = 1
        RETURN
      ENDIF

      IF(NTRACK_D.NE.1) THEN
        ICFAIL = 1
        RETURN
      ENDIF
      
      IF(NTRACK_D.EQ.1.AND.PZ_D(1).EQ.0.) THEN
        ICFAIL = 1
        RETURN
      ENDIF

      CALL ADCUNP('IC','CAL')
      CALL ADCUNP('RD','CMB')
      CALL ADCUNP('BV','CMB')
      CALL ADCUNP('BM','CAL')

      EB4 = 0.
      DO I=1,8
        DO J=1,2
          EB4 = EB4 + BMCAL(I,J,4)
        END DO
      END DO
      
      
      DO 666 ITRACK = 1,NTRACK_D

        TRKPAR(1) = RADIUS_D(ITRACK)
        TRKPAR(2) = XCEN_D(ITRACK)
        TRKPAR(3) = YCEN_D(ITRACK)

        CALL HFILL(IOFF+2,FLOAT(NPTSXY_D(ITRACK)),0.,1.)
        CALL HFILL(IOFF+3,FLOAT(NPTSZ_D(ITRACK)),0.,1.)

C-- See if any of the hits on the track are in the upper hemisphere.        
C-- IC is the pointer to the common UTC_CANDIDATE for this track

        IC = IPTRACK_D(ITRACK)

        NH = NHIT_CAND(IC)

        PHIW1 = PIBY2
        DO IHIT=1,NH
           
C-- IP is the pointer to the common UTC_XYZ for this hit

          IP = NP_CAND(IHIT,IC)

          ICELL = ICELL_XY(IP)
          ILAY  = ILAY_XY(IP)
          ISLAY = ISLAY_XY(IP)
          IF(AWPOS(ICELL,ILAY,ISLAY).GT.0..AND.
     &       AWPOS(ICELL,ILAY,ISLAY).LT.PI) THEN
            PHIW1 = AWPOS(ICELL,ILAY,ISLAY)
          ENDIF
        END DO

C-- Count the number of layers in the fit

        CALL VZERO(NHITPL,24)
        DO IHIT=1,NH
          IP = NP_CAND(IHIT,IC)
          ICELL = ICELL_XY(IP)
          ILAY  = ILAY_XY(IP)
          ISLAY = ISLAY_XY(IP)

          NLAYER = (ISLAY-1)/2*NLAY_MAX+ILAY   ! runs from 1 to 12

          DPHIW = PHIW1-AWPOS(ICELL,ILAY,ISLAY)
          IF(DPHIW.GT.PI) DPHIW=DPHIW-TWOPI
          IF(DPHIW.LT.-PI) DPHIW=DPHIW+TWOPI
          IF(ABS(DPHIW).GE.PIBY2) THEN    ! This hit is on the lower half
            NLAYER = NLAYER+12
          ENDIF

          NHITPL(NLAYER) = NHITPL(NLAYER)+1
        END DO
        
        NLAYFIT=0
        DO ILAYER=1,24
          IF(ILAYER.LE.12) THEN
            IF(NHITPL(ILAYER).GE.1.AND.USELAY(ILAYER)) NLAYFIT=NLAYFIT+1
          ELSE
            IF(NHITPL(ILAYER).GE.1.AND.USELAY(ILAYER-12))
     &                                            NLAYFIT=NLAYFIT+1
          ENDIF
        END DO

        CALL HFILL(IOFF+4,FLOAT(NLAYFIT),0.,1.)
        IF(NLAYFIT.LT.NXYPTS_CUT) THEN
          ICFAIL = 1
          RETURN
        ENDIF
        
        CALL HFILL(IOFF+5,CHARGE_D(ITRACK),0.,1.)

C-- Note that CHIXY_D is actually chisq/dof!!!!!!!

        IF(UTBFLD.NE.0.) THEN
          NDOF = NPTSXY_D(ITRACK)-3
        ELSE
          NDOF = NPTSXY_D(ITRACK)-2
        ENDIF
        
        CHIXY = CHIXY_D(ITRACK)*FLOAT(NDOF)
        
        IF(NDOF.GT.0.AND.CHIXY.GT.0.) THEN
          CHIPXY = PROB(CHIXY,NDOF)
          CHIDFXY = CHIXY_D(ITRACK)
        ELSE
          CHIPXY=-99.
          CHIDFXY=-99.
        ENDIF
        CALL HFILL(IOFF+6,CHIPXY,0.,1.)
        CALL HFILL(IOFF+7,CHIDFXY,0.,1.)

C-- Same for CHIZ_D

        NDOF = NPTSZ_D(ITRACK)-2
        CHIZ = CHIZ_D(ITRACK)*FLOAT(NDOF)
        IF(NDOF.GT.0.AND.CHIZ.GT.0.) THEN
          CHIPZ = PROB(CHIZ,NDOF)
          CHIDFZ = CHIZ_D(ITRACK)
        ELSE
          CHIPZ=-99.
          CHIDFZ=-99.
        ENDIF
        CALL HFILL(IOFF+8,CHIPZ,0.,1.)
        CALL HFILL(IOFF+9,CHIDFZ,0.,1.)

        PHI0 = ATAN2(COS2_D(ITRACK),COS1_D(ITRACK))
        IF(PHI0.LT.0.) PHI0 = PHI0+TWOPI

        CALL HFILL(IOFF+10,PHI0,DIST0_D(ITRACK),1.)
        CALL HFILL(IOFF+11,Z0_D(ITRACK),0.,1.)
        CALL HFILL(IOFF+12,Z0_D(ITRACK),DIST0_D(ITRACK),1.)

C-- Extrapolate the track to the I-counter

        IF(NPTSZ_D(ITRACK).GE.NHITZ_MIN) THEN
          CALL HFILL(IOFF+21,CHIDFZ,0.,1.)
          IF(CHIDFZ.LT.100.0) THEN
            DO I=1,2
              IF(I.EQ.2) THEN              ! Flip the track
                COS1_D(ITRACK) = -COS1_D(ITRACK)
                COS2_D(ITRACK) = -COS2_D(ITRACK)
                COS3_D(ITRACK) = -COS3_D(ITRACK)
              ENDIF
              CALL TRKPOS(ITRACK,RIC,0.,0.,XIC,YIC,ZIC,
     &                     COSIC(1),COSIC(2),COSIC(3))
              IF(XIC.NE.0..OR.YIC.NE.0.) THEN
                PHIIC = ATAN2(YIC,XIC)
                IF(PHIIC.LT.0.) PHIIC=PHIIC+TWOPI
                ISECIC = INT(3.0*PHIIC/PI) + 1      ! Ignore rotation
                EIC = ICCAL(ISECIC,1)
                CALL HFILL(IOFF+22,EIC,0.,1.)
                IF(EIC.GT.0.3) THEN
                  CALL HFILL(IOFF+23,ZIC,0.,1.)
                ENDIF
              ENDIF

C-- Extrapolate the track to the T-counter

              RTC=45.41                ! Radius of middle of T-counter
                                       ! Cannot use a parameter for flat_corr!

              CALL TRKPOS(ITRACK,RTC,0.,0.,XTC,YTC,ZTC,
     &                          COSTC(1),COSTC(2),COSTC(3))

              IF(XTC.NE.0..OR.YTC.NE.0.) THEN
                CALL FLAT_CORR(XTC,YTC,RTC)
                CALL TRKPOS(ITRACK,RTC,0.,0.,XTC,YTC,ZTC,
     &                           COSTC(1),COSTC(2),COSTC(3))
                PHITC = ATAN2(YTC,XTC)
                IF(PHITC.LT.0.) PHITC=PHITC+TWOPI
                ISECTC = INT(12.0*PHITC/PI) + 1
                ETC = RDCMB(1,1,ISECTC)
                CALL HFILL(IOFF+24,ETC,0.,1.)
                IF(ETC.GT.0.1) THEN
                  CALL HFILL(IOFF+25,ZTC,0.,1.)
                ENDIF
              ENDIF
              
C-- B4 extrapolation

              RB4=5.1                ! 'Radius' of B4 outline
              CALL TRKPOS(ITRACK,RB4,0.,0.,XB4,YB4,ZB4(I),
     &                      COSB4(1),COSB4(2),COSB4(3))

              IF(XB4.NE.0..OR.YB4.NE.0.) THEN
                CALL FLAT_CORR_B4(XB4,YB4,RB4)
                CALL TRKPOS(ITRACK,RB4,0.,0.,XB4,YB4,ZB4(I),
     &                      COSB4(1),COSB4(2),COSB4(3))
              ELSE
                ZB4(I) = 999.
              ENDIF
              
              IF(I.EQ.2) THEN              ! Flip the track back
                COS1_D(ITRACK) = -COS1_D(ITRACK)
                COS2_D(ITRACK) = -COS2_D(ITRACK)
                COS3_D(ITRACK) = -COS3_D(ITRACK)
              ENDIF

            END DO

            CALL HFILL(IOFF+26,EB4,0.,1.)
            IF(EB4.GT.0.3) THEN
              CALL HFILL(IOFF+27,MIN(ZB4(1),ZB4(2)),0.,1.)
            ELSE
              CALL HFILL(IOFF+28,MIN(ZB4(1),ZB4(2)),0.,1.)
            ENDIF

          ENDIF          ! Chiz/ndof
        ENDIF            ! nptsz_d
        
C-- Now histograms for hits
C-- IC is the pointer to the common UTC_CANDIDATE for this track

        IC = IPTRACK_D(ITRACK)

C-- Loop over the r-phi hits on the track

        NH = NHIT_CAND(IC)
        NHIT=0
        CALL VZERO(NHITPL,24)

        DO IHIT=1,NH

C-- IP is the pointer to the common UTC_XYZ for this hit
C-- IPTDC is the pointer to the TDC info (UTC_TDCDATA.CMN) for this hit

          IP = NP_CAND(IHIT,IC)
          IPTDC = IPTR_XY(IP)

          ICELL = ICELL_XY(IP)
          ILAY  = ILAY_XY(IP)
          ISLAY = ISLAY_XY(IP)
          NLAYER = (ISLAY-1)/2 * NLAY_MAX + ILAY
          JLAYER = NLAYER
          DPHIW = PHIW1-AWPOS(ICELL,ILAY,ISLAY)
          IF(DPHIW.GT.PI) DPHIW=DPHIW-TWOPI
          IF(DPHIW.LT.-PI) DPHIW=DPHIW+TWOPI
          IF(ABS(DPHIW).GE.PIBY2) THEN    ! This hit is on the lower half
            NLAYER = NLAYER+12
          ENDIF
          NHIT=NHIT+1
          NHITPL(NLAYER)=NHITPL(NLAYER)+1
          
          IF(NHIT.GT.MNHPT) THEN
            CALL KERROR(-1,10,'DPLOT','Too many rphi hits on a track')
            GO TO 666
          ENDIF
          IF(NHITPL(NLAYER).GT.MNHPL) THEN
            CALL KERROR(-1,10,'DPLOT','Too many rphi hits in a layer')
            GO TO 666
          ENDIF
            
          IWIRE(NHIT) = IWIRE_NUM(ICELL,ILAY,ISLAY)
          TIME(NHIT)  = LTDCC_U(IPTDC)
          LHIT(NHIT)  = .TRUE.
          PHW = UTC_WIRE_ADCC(ICELL,ILAY,ISLAY)

          IHITPL(NHITPL(NLAYER),NLAYER) = IHIT

          CALL HFILL(IOFF+101,RESID_CAND(IHIT,IC),FLOAT(JLAYER),1.)
          CALL HFILL(IOFF+102,RESID_S_CAND(IHIT,IC),FLOAT(JLAYER),1.)
          CALL HFILL(IOFF+103,RESID_CAND(IHIT,IC),PHW,1.)
          CALL HFILL(IOFF+104,PHW,FLOAT(JLAYER),1.)
          CALL HFILL(IOFF+105,RESID_CAND(IHIT,IC),
     &                                 FLOAT(NPTSXY_D(ITRACK)),1.)
          
          CALL CALCDCA(ICELL,ILAY,ISLAY,TRKPAR,XYPCA,DOCA,PHIVU,
     &                 PHIDCA)
          IF(PHIDCA.GT.0.) THEN
            DDRIFT = DRIFT_L(IP)
          ELSE
            DDRIFT = DRIFT_R(IP)
          ENDIF
          DOCA = SIGN(DOCA,-PHIDCA)
          DDRIFT = SIGN(DDRIFT,DOCA)
          RESID = DOCA-DDRIFT
          RESID2 = SIGN(RESID,RESID_CAND(IHIT,IC))

          CALL HFILL(IOFF+201,RESID,FLOAT(JLAYER),1.)
          CALL HFILL(IOFF+202,RESID2-RESID_CAND(IHIT,IC),0.,1.)

        END DO
C-- Loop over the Z hits on the track

        DO IHIT=1,NZCL(ITRACK)

C-- Loop over the strips in the cluster

          PHS = 0.
          DO IS = 1,NSIZCL(IHIT,ITRACK)
            IND = IPSZCL(IS,IHIT,ITRACK)
            ISTRIP = ISEC_NUM(IND)
            IFOIL  = IFOIL_NUM(IND)
            PHS = PHS + UTC_FOIL_ADCC(ISTRIP,IFOIL)
          END DO
            
          CALL HFILL(IOFF+111,RESZCL(IHIT,ITRACK),FLOAT(IFOIL),1.)
          CALL HFILL(IOFF+112,RESZCL(IHIT,ITRACK),PHS,1.)
          CALL HFILL(IOFF+113,PHS,FLOAT(IFOIL),1.)
          CALL HFILL(IOFF+114,RESZCL(IHIT,ITRACK),
     &                                FLOAT(NPTSZ_D(ITRACK)),1.)
          IF(NZCL(ITRACK).GE.NHITZ_MIN)
     &         CALL HFILL(IOFF+130+IFOIL,RESZCL(IHIT,ITRACK),
     &                                ZPOSCL(IHIT,ITRACK),1.)

C-- Figure out if the foil hit is on the upper or lower half of the
C-- detector

          IF(UTBFLD.EQ.0.) THEN
            IF( (PHI0.GT.PI.AND.PHICL(IHIT,ITRACK).LT.0.).OR.
     &          (PHI0.LT.PI.AND.PHICL(IHIT,ITRACK).GT.0.) ) THEN
              LFUP = .TRUE.
            ELSE
              LFUP = .FALSE.
            ENDIF
          ELSE
            IF( (PHI0.GT.PI.AND.
     &          -CHARGE_D(ITRACK)*PHICL(IHIT,ITRACK).LT.0.).OR.
     &          (PHI0.LT.PI.AND.
     &          -CHARGE_D(ITRACK)*PHICL(IHIT,ITRACK).GT.0.) ) THEN
              LFUP = .TRUE.
            ELSE
              LFUP = .FALSE.
            ENDIF
          ENDIF
             
C-- Loop over the wire hits in the adjacent layer

          IF(MOD(IFOIL,2).EQ.0) THEN
            NLAYER=2*IFOIL
          ELSE
            NLAYER=2*IFOIL-1
          ENDIF
          IF(.NOT.LFUP) NLAYER=NLAYER+12
          PHW=0.
          TMIN = 9999.
          ICELL_MIN = 0
          DO IW=1,NHITPL(NLAYER)
            IP = NP_CAND(IHITPL(IW,NLAYER),IC)
            ICELL = ICELL_XY(IP)
            ILAY  = ILAY_XY(IP)
            ISLAY = ISLAY_XY(IP)
            DPHIW = PHIW1-AWPOS(ICELL,ILAY,ISLAY)
            IF(DPHIW.GT.PI) DPHIW=DPHIW-TWOPI
            IF(DPHIW.LT.-PI) DPHIW=DPHIW+TWOPI
            IF(ABS(DPHIW).LT.PIBY2) THEN  ! This hit is on the upper half
              IF(LFUP) THEN               ! Foil hit is also on upper half
                PHW = PHW+UTC_WIRE_ADCC(ICELL,ILAY,ISLAY)
                IPTDC = IPTR_XY(IP)
                IF(LTDCC_U(IPTDC).LT.TMIN) THEN
                  TMIN = LTDCC_U(IPTDC)
                  ICELL_MIN = ICELL
                ENDIF
              ENDIF
            ELSE
              IF(.NOT.LFUP) THEN
                PHW = PHW+UTC_WIRE_ADCC(ICELL,ILAY,ISLAY)
                IPTDC = IPTR_XY(IP)
                IF(LTDCC_U(IPTDC).LT.TMIN) THEN
                  TMIN = LTDCC_U(IPTDC)
                  ICELL_MIN = ICELL
                ENDIF
              ENDIF
            ENDIF
          END DO

C-- Extrapolate the track to this foil and look at Z residual as a fn
C-- of strip. Note that for reasons of laziness we are excluding cases
C-- where there is no wire hit in the adjacent layer; in such cases
C-- we have to do a little extra computation to figure out which strip
C-- is intersected.

          IF(ICELL_MIN.GT.0) THEN
            CALL HFILL(IOFF+120+IFOIL,PHW,PHS,1.)
            SINL = COS3_D(ITRACK)
            COSL = SQRT(1.0-SINL*SINL)
            TANL = SINL/COSL
            IF(UTBFLD.EQ.0.) THEN
              IF(PHI0.GT.PI) THEN  ! track points down in r-phi
                ZFOIL = Z0_D(ITRACK) - PHICL(IHIT,ITRACK)*TANL
              ELSE
                ZFOIL = Z0_D(ITRACK) + PHICL(IHIT,ITRACK)*TANL
              ENDIF
            ELSE
              ZFOIL = Z0_D(ITRACK) + PHICL(IHIT,ITRACK)*TANL*
     &                            (-CHARGE_D(ITRACK))*RADIUS_D(ITRACK)
            ENDIF
               
            IF(ABS(ZFOIL).LE.0.5*ZLENL(ISLAY)) THEN
              DZMIN = 9999.
              DO IS = 1,NSIZCL(IHIT,ITRACK)
                IND = IPSZCL(IS,IHIT,ITRACK)
                ISTRIP = ISEC_NUM(IND)
                IF(FOIL_ZINT(ISTRIP,ICELL_MIN,IFOIL).GT.-60.) THEN
                  DZ = ABS(ZFOIL - FOIL_ZINT(ISTRIP,ICELL_MIN,IFOIL))
                  IF(DZ.LT.DZMIN) THEN
                    ISTRIP_MIN = ISTRIP
                    DZMIN = DZ
                  ENDIF
                ENDIF
              END DO
              IF(NZCL(ITRACK).GE.NHITZ_MIN) THEN
                  CALL HFILL(IOFF+140+IFOIL,RESZCL(IHIT,ITRACK),
     &                                     FLOAT(ISTRIP_MIN),1.) 
              ENDIF
            ENDIF
          ENDIF    
        END DO

666   CONTINUE                    ! End of loop over tracks

C-- Now fit the two branches separately

      COSMIC = .FALSE.
      CALL UTC_TRACK(.TRUE.,.TRUE.)
      
      IF(IFAIL_CODE(I_UT).NE.0) THEN
        ICFAIL = 1
        RETURN
      ENDIF

      CALL HFILL(IOFF2+1,FLOAT(NTRACK_D),0.,1.)
      
      IF(NTRACK_D.LT.1.OR.NTRACK_D.GT.2) THEN
        ICFAIL = 1
        RETURN
      ENDIF

      
C-- Now look at UTC-RSSC match

      IF(PARAM(1).NE.0.) THEN
        NTRACK_RS=0
        CALL RD_MAIN(2,0,'B')
        CALL HFILL(IOFF2+2,FLOAT(NTRACK_RS),0.,1.)
        DO I = 1,NTRACK_RS
          ITRACK = ITRACK_RS(I)
          IF(PZ_D(ITRACK).NE.0.) THEN
            CALL HFILL(IOFF2+3,FLOAT(NPTSXY_D(ITRACK)),0.,1.)
            CALL HFILL(IOFF2+4,FLOAT(NPTSZ_D(ITRACK)),0.,1.)
            IF(RS_X(I,1).LT.100..AND.RS_Y(I,1).LT.100..AND.
     &                             RS_Z(I,1).LT.100.) THEN
              CALL TRKPOS(ITRACK,0.,RS_X(I,1),RS_Y(I,1),
     &                  XSC1,YSC1,ZSC1,DXSC1,DYSC1,DZSC1)
              IF(XSC1.NE.0..OR.YSC1.NE.0.) THEN
                PHIUT1 = ATAN2(YSC1,XSC1)
                IF(PHIUT1.LT.0.) PHIUT1=PHIUT1+TWOPI
                PHISC1 = ATAN2(RS_Y(I,1),RS_X(I,1))
                IF(PHISC1.LT.0.) PHISC1=PHISC1+TWOPI
                DPHI1 = PHIUT1-PHISC1
                IF(DPHI1.GT.PI) DPHI1 = DPHI1-TWOPI
                IF(DPHI1.LT.-PI) DPHI1 = DPHI1+TWOPI
                DUTSC1 = SQRT( (RS_X(I,1)-XSC1)**2 +
     &                                  (RS_Y(I,1)-YSC1)**2 )
                DUTSC1 = SIGN(DUTSC1,DPHI1)
                ISEC1 = ISEC_PC(I,1)
                DZ1 = ZSC1 - RS_Z(I,1) ! - RSSC_ZOFF(ISEC1,1)
                CALL HFILL(IOFF2+301,DUTSC1,FLOAT(ISEC1),1.)
                CALL HFILL(IOFF2+303,DZ1,FLOAT(ISEC1),1.)
                CALL HFILL(IOFF2+305,DZ1,ZSC1,1.)
              ENDIF
            ENDIF

            IF(RS_X(I,2).LT.100..AND.RS_Y(I,2).LT.100..AND.
     &                             RS_Z(I,2).LT.100.) THEN
              CALL TRKPOS(ITRACK,0.,RS_X(I,2),RS_Y(I,2),
     &                  XSC2,YSC2,ZSC2,DXSC2,DYSC2,DZSC2)
             
              IF(XSC2.NE.0..OR.YSC2.NE.0.) THEN
                PHIUT2 = ATAN2(YSC2,XSC2)
                IF(PHIUT2.LT.0.) PHIUT2=PHIUT2+TWOPI
                PHISC2 = ATAN2(RS_Y(I,2),RS_X(I,2))
                IF(PHISC2.LT.0.) PHISC2=PHISC2+TWOPI
                DPHI2 = PHIUT2-PHISC2
                IF(DPHI2.GT.PI) DPHI2 = DPHI2-TWOPI
                IF(DPHI2.LT.-PI) DPHI2 = DPHI2+TWOPI
                DUTSC2 = SQRT( (RS_X(I,2)-XSC2)**2 +
     &                                  (RS_Y(I,2)-YSC2)**2 )
                DUTSC2 = SIGN(DUTSC2,DPHI2)
                ISEC2 = ISEC_PC(I,2)
                DZ2 = ZSC2 - RS_Z(I,2) ! - RSSC_ZOFF(ISEC2,2)
                CALL HFILL(IOFF2+302,DUTSC2,FLOAT(ISEC2),1.)
                CALL HFILL(IOFF2+304,DZ2,FLOAT(ISEC2),1.)
                CALL HFILL(IOFF2+306,DZ2,ZSC2,1.)
              ENDIF
            ENDIF

          ENDIF
        END DO
      ENDIF
      
C-- Now look at correlations between the two legs of the cosmic ray track

      IF(NTRACK_D.EQ.2.AND.PZ_D(1).NE.0..AND.PZ_D(2).NE.0.) THEN
         
        SINL1 = COS3_D(1)
        COSL1 = SQRT(1.0-SINL1*SINL1)
        TANL1 = SINL1/COSL1

        SINL2 = COS3_D(2)
        COSL2 = SQRT(1.0-SINL2*SINL2)
        TANL2 = SINL2/COSL2

        PHI1 = ATAN2(COS2_D(1),COS1_D(1))
        IF(PHI1.LT.0.) PHI1 = PHI1+TWOPI

        PHI2 = ATAN2(COS2_D(2),COS1_D(2))
        IF(PHI2.LT.0.) PHI2 = PHI2+TWOPI

        Z01 = Z0_D(1)
        Z02 = Z0_D(2)

        D01 = DIST0_D(1)
        D02 = DIST0_D(2)

        CALL TRKPOS(1,RFOIL(6),0.,0.,XFOIL,YFOIL,ZFOIL,
     &                                DXFOIL,DYFOIL,DZFOIL)
        Z1 = ZFOIL
          
        CALL TRKPOS(2,RFOIL(6),0.,0.,XFOIL,YFOIL,ZFOIL,
     &                                DXFOIL,DYFOIL,DZFOIL)
        Z2 = ZFOIL

C-- Define track 1 to be the one with positive tan(lambda)
C-- If both tracks are positive, then this is a random selection

        IF(SINL1.LT.0.) THEN
          DUMMY = TANL1
          TANL1 = TANL2
          TANL2 = DUMMY

          DUMMY = PHI1
          PHI1 = PHI2
          PHI2 = DUMMY

          DUMMY = Z01
          Z01 = Z02
          Z02 = DUMMY

          DUMMY = D01
          D01 = D02
          D02 = DUMMY

          DUMMY = Z1
          Z1 = Z2
          Z2 = DUMMY
        ENDIF

        CALL HFILL(IOFF2+211,TANL1-ABS(TANL2),TANL1,1.)
        CALL HFILL(IOFF2+212,Z01-Z02,Z01,1.)

        IF(UTBFLD.EQ.0.) THEN
          DELTAD = D01+D02
        ELSE
          DELTAD = D01-D02
        ENDIF
        CALL HFILL(IOFF2+214,DELTAD,D01,1.)

        DELTAP = PHI1-PHI2
        IF(DELTAP.GT.0.) THEN
          DELTAP = DELTAP - PI
        ELSE
          DELTAP = DELTAP + PI
        ENDIF

        CALL HFILL(IOFF2+213,DELTAP,PHI1,1.)
        CALL HFILL(IOFF2+215,Z01-Z02,TANL1,1.)

C-- Look for a twist in the chamber
C-- Select back-to-back cosmics

        IF( (SINL1.GT.0..AND.SINL2.LT.0.).OR.
     &      (SINL1.LT.0..AND.SINL2.GT.0.) ) THEN
          IF(TANL1+TANL2.LT.0.07.AND.ABS(Z01-Z02).LT.1.0) THEN
            CALL HFILL(IOFF2+216,DELTAP,Z1-Z2,1.)
            CALL HFILL(IOFF2+217,DELTAD,Z1-Z2,1.)
          ENDIF
        ENDIF

C-- Check for a coherent miscalibration of the foil stereo angles

        IF( (Z1.GT.0.0.AND.Z2.GT.0.0) .OR.
     &      (Z1.LT.0.0.AND.Z2.LT.0.0) ) THEN
          CALL HFILL(IOFF2+218,TANL1-ABS(TANL2),TANL1,1.)
        ENDIF

        IF( (Z1.GT.0.0.AND.Z2.LT.0.0) .OR.
     &      (Z1.LT.0.0.AND.Z2.GT.0.0) ) THEN
          CALL HFILL(IOFF2+219,TANL1-ABS(TANL2),TANL1,1.)
        ENDIF

      ENDIF
      
C-- Fit the two branches separately again. This time leave the
C-- middle superlayer out of the fit

      DO ILAY=5,8
        USELAY(ILAY) = .FALSE.
      END DO

      CALL UTC_TRACK(.TRUE.,.TRUE.)

      IF(NTRACK_D.LT.1.OR.NTRACK_D.GT.2) THEN
        ICFAIL = 1
        RETURN
      ENDIF

C-- Look at signed DOCA residuals as a function of Z

      DO ITRACK=1,NTRACK_D
        IF(PZ_D(ITRACK).NE.0.0.AND.NPTSZ_D(ITRACK).GE.4) THEN
          IC = IPTRACK_D(ITRACK)
          NH = NHIT_CAND(IC)
          PHI0 = ATAN2(COS2_D(ITRACK),COS1_D(ITRACK))
          IF(PHI0.LT.0.) PHI0 = PHI0+TWOPI

          DO IHIT=1,NH
            IP = NP_CAND(IHIT,IC)
            ICELL = ICELL_XY(IP)
            ILAY  = ILAY_XY(IP)
            ISLAY = ISLAY_XY(IP)
            NLAYER = (ISLAY-1)/2 * NLAY_MAX + ILAY

            IF(NLAYER.GE.5.AND.NLAYER.LE.8) THEN
              RESID = RESID_CAND(IHIT,IC)
              ILR = ILR_CAND(IHIT,IC)
              IF(ILR.EQ.2) RESID = -RESID
              CALL TRKPOS(ITRACK,0.,XWPOS(ICELL,ILAY,ICELL),
     &          YWPOS(ICELL,ILAY,ISLAY),X,Y,Z,DX,DY,DZ)
              IF(X.NE.0..OR.Y.NE.0.) THEN
                CALL HFILL(IOFF2+101,Z,RESID,1.)
                CALL HFILL(IOFF2+102+NLAYER-5,Z,RESID,1.)
                IF(ABS(Z).LT.25.0) THEN
                  IZBIN = INT((Z+25.)/5.) + 1
                  CALL HFILL(IOFF2+110+IZBIN-1,PHI0,RESID,1.)
                ENDIF
              ENDIF
            ENDIF
          END DO
        ENDIF
      END DO
      
      DO ILAY=5,8
        USELAY(ILAY) = .TRUE.
      END DO
      
C-- Fit the two RSSC branches to a single track in Z. The math below
C-- assumes that the B field is 0.

      IF(NTRACK_RS.EQ.2.AND.UTBFLD.EQ.0.) then
        ITRK1 = ITRACK_RS(1)
        ITRK2 = ITRACK_RS(2)
        IF(PZ_D(ITRK1).EQ.0..OR.PZ_D(ITRK2).EQ.0.) GO TO 7734
        
        SINL1 = COS3_D(ITRK1)
        COSL1 = SQRT(1.0-SINL1*SINL1)
        TANL1 = SINL1/COSL1

        SINL2 = COS3_D(ITRK2)
        COSL2 = SQRT(1.0-SINL2*SINL2)
        TANL2 = SINL2/COSL2

        PHI1 = ATAN2(COS2_D(ITRK1),COS1_D(ITRK1))
        IF(PHI1.LT.0.) PHI1 = PHI1+TWOPI

        PHI2 = ATAN2(COS2_D(ITRK2),COS1_D(ITRK2))
        IF(PHI2.LT.0.) PHI2 = PHI2+TWOPI

        DELTAP = PHI1-PHI2
        IF(DELTAP.GT.0.) THEN
          DELTAP = DELTAP - PI
        ELSE
          DELTAP = DELTAP + PI
        ENDIF

        Z01 = Z0_D(ITRK1)
        Z02 = Z0_D(ITRK2)

        D01 = DIST0_D(ITRK1)
        D02 = DIST0_D(ITRK2)
        DELTAD = D01+D02

C-- Select shallow dip angle

        TLMAX = MAX(ABS(TANL1),ABS(TANL2))
c        CALL HFILL(IOFF2+310,TLMAX,0.,1.)
        IF(TLMAX.GT.0.3) GO TO 7734
        
C-- Select back-to-back cosmics

        IF(ABS(TANL1+TANL2).GT.0.03) GO TO 7734

        IF(ABS(Z01-Z02).GT.0.6) GO TO 7734

        IF(ABS(DELTAP).GT.0.005) GO TO 7734

        IF(ABS(DELTAD).GT.0.07) GO TO 7734
        
C-- Require 21 layers in the range stack on both branches

        IMOD = ITRK(NHTRK(1),1)
        ISEC = 1+(IMOD-1)/21
        ILAY = IMOD-21*(ISEC-1)
        IF(ILAY.NE.21) GO TO 7734

        IMOD = ITRK(NHTRK(2),2)
        ISEC = 1+(IMOD-1)/21
        ILAY = IMOD-21*(ISEC-1)
        IF(ILAY.NE.21) GO TO 7734
        
C-- Look at BV energy
        
        RBV=89.59                ! Radius of back of layer 21
        CALL TRKPOS(ITRK1,RBV,0.,0.,XBV,YBV,ZBV,
     &                          COSBV(1),COSBV(2),COSBV(3))
        IF(XBV.NE.0..OR.YBV.NE.0.) THEN
          CALL FLAT_CORR(XBV,YBV,RBV)
          CALL TRKPOS(ITRK1,RBV,0.,0.,XBV,YBV,ZBV,
     &                           COSBV(1),COSBV(2),COSBV(3))
          PHIBV = ATAN2(YBV,XBV)
          IF(PHIBV.LT.0.) PHIBV=PHIBV+TWOPI
          ISECBV = INT(12.0*PHIBV/PI)+1
          ISECBV = 2*ISECBV
          EBV = 0.
          EBVLAST = 0.
          DO ISEC=ISECBV-3,ISECBV+2
            JSEC = ISEC
            IF(ISEC.LT.0) JSEC=ISEC+48
            IF(ISEC.GT.48) JSEC=ISEC-48
            DO ILAY=1,4
              EBV = EBV+BVCMB(1,ILAY,JSEC)
            END DO
            EBVLAST = EBVLAST+BVCMB(1,4,JSEC)
          END DO
          IF(EBV.LT.50.0.OR.EBV.GT.100.0) GO TO 7734
          IF(EBVLAST.LT.12.0.OR.EBVLAST.GT.30.0) GO TO 7734
        ENDIF
        
        RBV=89.59                ! Radius of back of layer 21
        CALL TRKPOS(ITRK2,RBV,0.,0.,XBV,YBV,ZBV,
     &                          COSBV(1),COSBV(2),COSBV(3))
        IF(XBV.NE.0..OR.YBV.NE.0.) THEN
          CALL FLAT_CORR(XBV,YBV,RBV)
          CALL TRKPOS(ITRK2,RBV,0.,0.,XBV,YBV,ZBV,
     &                           COSBV(1),COSBV(2),COSBV(3))
          PHIBV = ATAN2(YBV,XBV)
          IF(PHIBV.LT.0.) PHIBV=PHIBV+TWOPI
          ISECBV = INT(12.0*PHIBV/PI)+1
          ISECBV = 2*ISECBV
          EBV = 0.
          EBVLAST = 0.
          DO ISEC=ISECBV-3,ISECBV+2
            JSEC = ISEC
            IF(ISEC.LT.0) JSEC=ISEC+48
            IF(ISEC.GT.48) JSEC=ISEC-48
            DO ILAY=1,4
              EBV = EBV+BVCMB(1,ILAY,JSEC)
            END DO
            EBVLAST = EBVLAST+BVCMB(1,4,JSEC)
          END DO
          IF(EBV.LT.50.0.OR.EBV.GT.100.0) GO TO 7734
          IF(EBVLAST.LT.12.0.OR.EBVLAST.GT.30.0) GO TO 7734
        ENDIF
        
        DO I=1,4
          IF(I.LE.2) THEN
            IF(RS_X(1,I).GT.100..OR.RS_Y(1,I).GT.100..OR.
     &                             RS_Z(1,I).GT.100.) GO TO 7734
            RSC(I) = RS_X(1,I)**2 + RS_Y(1,I)**2
            PATH = RSC(I) - DIST0_D(ITRK1)**2
            IF(PATH.LT.0.) GO TO 7734
            PATH = SQRT(PATH)
            DPHI(I) = PATH
            ISECT(I) = ISEC_PC(1,I)
            ZED(I) = RS_Z(1,I)! + RSSC_ZOFF(ISECT(I),I)
          ELSE
            IF(RS_X(2,I-2).GT.100..OR.RS_Y(2,I-2).GT.100..OR.
     &                             RS_Z(2,I-2).GT.100.) GO TO 7734
            RSC(I) = RS_X(2,I-2)**2 + RS_Y(2,I-2)**2
            PATH = RSC(I) - DIST0_D(ITRK2)**2
            IF(PATH.LT.0.) GO TO 7734
            PATH = SQRT(PATH)
            DPHI(I) = -PATH
            ISECT(I) = ISEC_PC(2,I-2)
            ZED(I) = RS_Z(2,I-2)! + RSSC_ZOFF(ISECT(I),I-2)
          ENDIF
        END DO

        CALL COSMIC_LINFIT(4,LFIT,DPHI,ZED,1.5,PARAMS,CHIZ,NDOF,COV,
     &            DUMMY,IDUMMY,ZRESID,XTRAER)
        IF(NDOF.LT.2) GO TO 7734


        CHIPZ = PROB(CHIZ,NDOF)

        IF(CHIPZ.LT.0.01) GO TO 7734
        
        TANL = PARAMS(1)
        PHI0 = ATAN2(COS2_D(ITRK1),COS1_D(ITRK1))
        IF(PHI0.LT.0.) PHI0 = PHI0 + TWOPI
        SINP = SIN(PHI0)
        COSP = COS(PHI0)
        DENOM = SQRT(1.0+TANL*TANL)
        SINTH = 1.0/DENOM
        COSTH = TANL/DENOM
        COS1_D(ITRK1) = SINTH*COSP
        COS2_D(ITRK1) = SINTH*SINP
        COS3_D(ITRK1) = COSTH
        Z0_D(ITRK1) = PARAMS(2)

        COS1_D(ITRK2) = -COS1_D(ITRK1)
        COS2_D(ITRK2) = -COS2_D(ITRK1)
        COS3_D(ITRK2) = -COS3_D(ITRK1)
        Z0_D(ITRK2) = Z0_D(ITRK1)

        DO I=1,2
          ITRACK = ITRACK_RS(I)
          CALL TRKPOS(ITRACK,RIC,0.,0.,XIC,YIC,ZIC,COSIC(1),COSIC(2),
     &              COSIC(3))
          IF(XIC.NE.0..OR.YIC.NE.0.) THEN
            PHIIC = ATAN2(YIC,XIC)
            IF(PHIIC.LT.0.) PHIIC=PHIIC+TWOPI
            ISECIC = INT(3.0*PHIIC/PI) + 1
            EIC = ICCAL(ISECIC,1)
          ENDIF

          RTC=45.41                ! Radius of middle of T-counter
          CALL TRKPOS(ITRACK,RTC,0.,0.,XTC,YTC,ZTC,COSTC(1),COSTC(2),
     &              COSTC(3))
          IF(XTC.NE.0..OR.YTC.NE.0.) THEN
            CALL FLAT_CORR(XTC,YTC,RTC)
            CALL TRKPOS(ITRACK,RTC,0.,0.,XTC,YTC,ZTC,COSTC(1),COSTC(2),
     &                COSTC(3))
            PHITC = ATAN2(YTC,XTC)
            IF(PHITC.LT.0.) PHITC=PHITC+TWOPI
            ISECTC = INT(12.0*PHITC/PI) + 1
            ETC = RDCMB(1,1,ISECTC)
          ENDIF
        END DO

      ENDIF
7734  CONTINUE
      
      RETURN

      END


      SUBROUTINE FLAT_CORR_B4(X,Y,RSEG)

C      correct for flat surfaces

      PI = 3.141592653 

      PHI = PI/2.                               
      IF( X .NE. 0. ) PHI = ATAN(ABS(Y/X))
      PHI = AMOD(PHI,.7854)
      PHI = ABS(PHI-.3927)
      RSEG = RSEG/COS(PHI)

      RETURN
      END

