+PATCH,$CORR.
+DECK,CORR1.                     Updates version:  2.01/00  to   2.01/14
+REP,*TITLE*,TITLE,1.
LOCATER   2.01/14  07/08/97  00.15.01
+REP,$VERSION,V2_01,1-2.
*CMZU:  2.01/14 07/08/97  00.15.01  by  Rafik Ouared
*CMZU:  2.01/13 29/03/97  21.56.16  by  Rafik Ouared
*CMZU:  2.01/12 03/03/97  16.15.31  by  Rafik Ouared
*CMZU:  2.01/11 14/02/97  09.25.42  by  Rafik Ouared
*CMZU:  2.01/10 26/01/97  14.16.18  by  Rafik Ouared
*CMZU:  2.01/09 26/11/96  20.08.33  by  Rafik Ouared
*CMZU:  2.01/08 27/06/96  18.10.47  by  F.-H.Heinsius
*CMZ :  2.01/07 17/06/96  18.09.04  by  F.-H.Heinsius
*CMZU:  2.01/06 03/06/96  11.01.29  by  Unknown
*CMZU:  2.01/05 02/04/96  19.54.02  by  F.-H.Heinsius
*CMZ :  2.01/04 20/03/96  18.03.38  by  F.-H.Heinsius
*CMZ :  2.01/03 15/03/96  18.31.02  by  F.-H.Heinsius
*CMZ :  2.01/02 14/03/96  16.48.58  by  F.-H.Heinsius
*CMZ :  2.01/01 13/03/96  15.39.42  by  F.-H.Heinsius
*CMZU:  2.01/00 21/12/95  18.20.11  by  Christian Voelcker
*-- Author :    Christian Voelcker   21/12/95
*________________________________________________________________________
*
*::>          VERSION  2.01/14   07/08/97 00.15.01
*
*             TCTRAK modified to avoid conflict between TCTK and TCTR
*             which is seen in GTSOFF.
*             TCCALB, TCFFUZ,TCFCIRD modified to avoid dependence of
*             the reconstruction upon the previous event.
*________________________________________________________________________
*
*::>          VERSION  2.01/13   29/03/97 21.56.15
*             Removed redudant code in TCINIB
*             Fix MZLIFT garbage collection errors in TCVER3,TCVRHX,TCVRTX
*________________________________________________________________________
*
*::>          VERSION  2.01/12   03/03/97 16.15.31
*
*             New routine TCINIB to read in SV calibration. (M.Doser)
*             TCFORM modified to allow transportability of RSVX bank
*             through different platforms.
*             plus minor modifications in TSTVTX,RSRVTX.
*________________________________________________________________________
*
*::>          VERSION  2.01/11   14/02/97 09.25.42
*
*             new fix on TCVRHX to prevent zero division (M.Lakata)
*________________________________________________________________________
*
*::>          VERSION  2.01/10   26/01/97 14.16.18
*             TCVER3 and TCVRHX modified by M.Lakata to prevent
*             handling extremely small and large numbers in these routines.
*             Besides, INSTALL and CJINIT have been updated either.
*________________________________________________________________________
*
*::>          VERSION  2.01/09   26/11/96 20.08.33
*             Include the SVX calibration constants, both
*             for global as for detector specific corrections
*             max number of hits in TCCIRC set to 50
*             SPP elements in TCVRHX are less than 1.D+10
*
*________________________________________________________________________
*
*::>          VERSION  2.01/08   27/06/96 18.10.45
*
*             Do certain garfield corections not for MC. Smearing is
*             done in position not in time (were drived for that).
*             Change mean crate offset from -5.5ns to +2.5 ns.
*             See: TJDCGT TJTIME TSMINI TSMUSH
*             Affects mainly error in 1/pxy.
*
*::>          VERSION  2.01/07   17/06/96 18.09.04
*
*             vertex fit: change GAXTC Parameter (TCVRHX) and reduce CPU time
*             by calling TCVERT for simple cases in TCVER3.
*
*::>          VERSION  2.01/06   03/06/96 11.01.28
*               update of installation script for solaris
*
*::>          VERSION  2.01/05   02/04/96 19.54.02
*
*             Complete modified TCVER3 vertex fit and updated
*             assignments of errors (M.Lakata).
*             Old versions have wrong errors.
*             INSTALL: remove default path from compile script,
*             i.e. add P option in SET XA... (FHH).
*
* List of modified decks
* Patch $KUMACS:  INSTALL
* Patch TCCOMMON: C_MACRO
* Patch TC_MAIN:  CJINIT
* Patch TC_RAWS:  TPWPOS TSRPHI
* Patch TC_RADON1:RADON1
* Patch TC_CIRC:  TXTGET
* Patch TC_HELX:  TPCNCT
* Patch TC_VERT:  TCVER3 TCVRHX TCVPCI TCVPCR
* Patch TC_UTIL:  TCPUSH TCPOP
*
*::>          VERSION  2.01/04   20/03/96 18.03.38
*
*             INSTALL script from cboff (M.Lakata)
*             TPWPOS: positive errors for PWC (M.Lakata)
*
*::>          VERSION  2.01/03   15/03/96 18.31.02
*
*             TCHXLD: Correct error in z=3.5cm for SVX.
*             INSTALL: New option INC (for include statements)
*
*::>          VERSION  2.01/02   14/03/96 16.48.57
*
*             New arrays added up in TCLIFT: MTSV1,MTSV2
*             TCFORM,TSVPOS modified to fix up bugs in
*             zebra bank allocations and to add a forth bank
*             in tsvpos. (R.Ouared)
*             TCFILL: nr_tracks_founds not filled if no hits
*             in event. (F.H.Heinsius)
*
*::>          VERSION  2.01/01   13/03/96 15.39.41
*
*             Updates in TCVER3 vertex fit from M.Lakata:
*             TXTGET: use error instead of error**2 for PWC/SVX
*             TCVRHZ: use IQ( +4) instead of Q()
*
*             Fix format statements in TCVRHX.
+REP,$KUMACS,HELP,1-32.
*CMZU:  2.01/04 18/03/96  16.05.02  by  F.-H.Heinsius
*--  Obsolete
*CMZ :  2.01/03 15/03/96  13.48.57  by  F.-H.Heinsius
*CMZU:  1.99/11 06/07/95  23.59.28  by  Curtis A. Meyer
*CMZU:  1.96/03 10/04/95  15.53.43  by  Curtis A. Meyer
*CMZU:  1.94/02 02/02/95  08.47.08  by  Curtis A. Meyer
*-- Author :    Curtis A. Meyer  02.Nov.1994
+ADD,$KUMACS,INSTALL,0.
*CMZU:  2.01/10 22/01/97  16.31.02  by  M.T. Lakata
*CMZU:  2.01/06 15/05/96  13.50.00  by  F.-H.Heinsius
*CMZU:  2.01/05 02/04/96  18.46.59  by  F.-H.Heinsius
*CMZ :  2.01/04 18/03/96  17.11.11  by  F.-H.Heinsius
*CMZ :  2.01/03 15/03/96  12.50.00  by  F.-H.Heinsius
+ADD,$KUMACS,INSTALL,29.
*   16/03/96 Mark Lakata
*   Cleaned up install so that it works like cboff
*
*
+REP,$KUMACS,INSTALL,35-71.
*  goto HELP for instructions
*
   if [src]='X' Then
      goto help
   Endif

   If [src]='?' Then
      goto help
   Endif
*
*  SRC    Source code directory
*
   if $os .eq. unix then
     if ( .not. ($FEXIST([SRC]) =1)) then
       mess Directory [SRC] does not exist.
       exitm
     endif
   endif

   message Putting source code files in [SRC]
*
*  set defaults
*
   BUILDLEVEL = LIB
   OBJECTSTAT = FALSE
   CCODEFLAG  = TRUE
   DEBUGFLAG  = FALSE
   GNUFLAG    = FALSE
   INCFLAG    = FALSE
   INTERNAL   = FALSE
   locversion    = ''
   cboffversion    = ''

   Select  .
   Select CCODE
+REP,$KUMACS,INSTALL,77-133.
     Case $Upper([%I]) In
       (DEBUG)       DEBUGFLAG = TRUE
                     message ' compiling with debug flag (-g -C). '
       (ALL)         BUILDLEVEL = LIB
                     OBJECTSTAT = KEEP
       (LIB)         BUILDLEVEL = LIB
       (CAR)         car locater.car
                     exitm
       (SOURCE)      BUILDLEVEL = SOURCE
       (OBJECT)      BUILDLEVEL = OBJECT
                     OBJECTSTAT = KEEP
       (CALIB)       Select JDCALIBR
                     message ' Calibration code has been selected. '
       (INC)         INCFLAG   = TRUE
                     message '* Creating #include files'
       (GNU)         GNUFLAG   = TRUE
                     message ' Gnu Compilers being used. '
       (DBGPRNT)     Select DEBFALSE
                     Select DEBPRINT
                     Select DEBUGPAT
                     Select DEBUGFIT
                     Select DEBHELIX
                     Select DEBVERTX
                     message ' Debug print lines activated '
       (DBGPRNTX)    Select DEBFALSE
                     Select DDDEBUG
                     Select DEBPRINT
                     Select DEBUGRAW
                     Select DEBUGPAT
                     Select DEBUGFIT
                     Select DEBHELIX
                     Select DEBVERTX
                     message ' Extreme Debug print lines activated '
       (BXIMPROVE)   Select BXIMPROVE
                     message ' BXIMPROVE selected '
       (BXUSER)      Select BXUSER
                     message ' BXUSER selected '
       (NOCCODE)     CCODEFLAG = FALSE
                     select -CCODE
                     message ' Deselecting C code '
       (HISTOG)      Select TCHBOOK
                     message ' Fast Fuzzy Histograms Selected '
       (INTERNAL)    Select INTERNAL
                     INTERNAL = TRUE
                     message ' Internal Locater commons selected '
       (LOCVERSION:*)   locversion = $SUBSTRING([%I],12,100)
                     message '* Extracting Locater version '//[locversion]
       (CBOFFVERSION:*) cboffversion = $SUBSTRING([%I],14,100)
                     message '* Extracting CBOFF version '//[cboffversion]
       (*)     message '****************************************'
               message Error: Option $Upper([%I]) not understood.
               message '****************************************'
               FAIL = TRUE
      Endcase
   Enddo

   If [FAIL] = TRUE then
     goto help
   Endif
   if ([locversion] .ne. '' .and. [cboffversion] = '' ) .or. ([locversion] = '' .and. [cboffversion] .ne. '') then
     mess 'You must specific both LOCVERSION:x.xx/xx and CBOFFVERSION:y.yy/yy'
     exitm
   endif

   Message * Build Level = [BUILDLEVEL]

   Case $OS In
     (UNIX) cboffcmz = $CBCMZ/cboff.cmz
     (VMS)  cboffcmz = CB$CMZ:cboff.cmz
   Endcase

   if ($fexist([cboffcmz]) = 1) then
       file [cboffcmz] -R
   else
       mess can not find [cboffcmz], check CBCMZ env var
       exitm
   endif

+ADD,$KUMACS,INSTALL,145.
     (SOLARIS) Select SUN
     (LINUX) Select SUN
+REP,$KUMACS,INSTALL,155-179.
   cext = c
   ext  = f
   Case $MACHINE In
      (IBM)    ext = fortran
      (VAX)    ext = for
      (ALPHA)  Case $OS In
                 (VMS)  ext = for
               Endcase
   Endcase
*
   Set FORTRAN -LAN
   if [OBJECTSTAT] = KEEP Then
      XACL = -XACLP
   else
      XACL = -XADCLP
   endif

   IF [INCFLAG] = TRUE Then
     Case $OS In
       (UNIX) set '#include "$SEQ.inc"'   -INC F77
              FileCase Keep
              ext = F
       (VMS)  set '      INCLUDE 'CB$INC:$SEQ.INC''   -INC F77
     Endcase
   Else
     set . -inc
   EndIf
*
* Set up filenames for the extraction of the fortran code.
*
   Case $OS In
     (UNIX) Set [SRC]/mklocater.sh    [XACL]
            Set [SRC]/*.[ext]         -F F77
            Set [SRC]/*.[cext]        -F C
            Set liblocater.a          -L
     (VMS)  Set [SRC]mklocater.com    [XACL]
            Set [SRC]*.[ext]          -F F77
            Set [SRC]*.[cext]         -F C
            Set LOCATER.OLB           -L
+REP,$KUMACS,INSTALL,188.
         (HPUX,SGI,SUN,DECS,IBMRT,SOLARIS,LINUX) DEBUGOPT=' -g -C'
+REP,$KUMACS,INSTALL,210-276.
   Case $MACHINE In
     (HPUX) Set 'f77 -c'//[DEBUGOPT]//' +ppu $compfile' -C f77
          Set 'cc -c -O -w '//[DEBUGCCC]//'        $compfile' -C C
     (LINUX)       Set 'g77 -c'//[DEBUGOPT]//' -w $compfile' -C f77
                   Set 'gcc -c'//[DEBUGCCC]//' -w $compfile' -C C
     (SUN,SOLARIS) Set 'f77 -c'//[DEBUGOPT]//' -w $compfile' -C f77
                   Set 'gcc -c'//[DEBUGCCC]//' -w $compfile' -C C
     (IBMRT) Set 'xlf -c'//[DEBUGOPT]//' -qextname $compfile' -C f77
           Set 'cc -c -w '//[DEBUGCCC]//'        $compfile' -C C
     (SGI)   Set 'f77 -c'//[DEBUGOPT]//' $compfile G 3' -C f77
           Set 'cc -c -w '//[DEBUGCCC]//'        $compfile' -C C
     (DECS)  Set 'f77 -c -w -Nl99 '//[DEBUGOPT]//' $compfile' -C f77
           Set 'cc -c -w '//[DEBUGCCC]//'        $compfile' -C C
     (VAX)   Set 'FOR/OBJ=$compfile.OBJ'//[DEBUGOPT]//' $compfile' -C f77
           Set 'CC/OBJ=$compfile.obj'//[DEBUGCCC]//'  $compfile' -C C
*
     (ALPHA) Case $OS In
      (VMS) Set 'FOR/SEPAR/OBJ=$compfile.OBJ'//[DEBUGOPT]//' $compfile' -C f77
          Set 'CC/OBJ=$compfile.obj'//[DEBUGCCC]//' $compfile' -C C
      (UNIX) Set 'f77 -c -O -w -Nl99 '//[DEBUGOPT]//' $compfile' -C f77
           Set 'cc -c -O -w '//[DEBUGCCC]//' $compfile' -C C
     Endcase
*
   Endcase
*
   IF [GNUFLAG] = TRUE THEN
     Set '`which g77` -c'//[DEBUGOPT]//' -w $compfile' -C f77
     Set '`which gcc` -c'//[DEBUGCCC]//' -w $compfile' -C C
   Endif
*
   IF [INTERNAL] = TRUE THEN
     if [locversion] = '' then
       Sequence //locater/tccommon
     else
       Sequence -V [locversion] //locater/tccommon
     endif
   Endif
   if [cboffversion] = '' then
     Sequence //cboff/commcb
   else
     Sequence -V [cboffversion] //cboff/commcb
   endif
*
* DO FORTRAN
*
   set F77 -LAN
   Cd //locater
   Case [BUILDLEVEL] In
     (NONE)
       message * no code extracted
     (SOURCE)
       if [locversion] = '' then
         Cxtract *
       else
         Cxtract -V [locversion] *
       endif
       Message * fortran code extracted
     (OBJECT)
       if [locversion] = '' then
         Ccompile *
       else
         Ccompile -V [locversion] *
       endif
       Message * fortran objects created
     (LIB)
       if [locversion] = '' then
         Clib *
       else
         Clib -V [locversion] *
       endif
       Message * fortran code added to library
   EndCase
*
   IF [CCODEFLAG] = TRUE Then
*
     Set C -LAN
*
     Cd //locater/tc_radon2
     Case [BUILDLEVEL] In
          (NONE)
            message 'no code extracted'
          (SOURCE)
            Cxtract *
            Message * c extracted
          (OBJECT)
            CCompile *
            Message * c object files compiled
          (LIB)
            CLib *
            Message * c code added to library
     EndCase
*
   Endif

   if [incflag] = true then
     mess ********** ! ! WARNING ! !  ***************************************
     mess The *.inc files may be duplicated in CBINC, and there is
     mess no guarantee that these include files match them.
     mess It is your responsibility to make sure that the include files
     mess are the same, otherwise you will have serious common block
     mess misalignment.
     mess The suggested method is to ONLY extract the code, then delete
     mess all these locally *.inc files, and then
     mess compile the code outside of CMZ using the -I$(CBINC) flag to
     mess verify that all the cb software uses the same include files found
     mess in CBINC.
     mess *******************************************************************
   endif
*
END:
*
 EXITM

HELP:
   mess 'Command template:'
   mess 'UNIX: cmz -install locater <files_directory> [ <option> ... ]'
   mess 'VMS:  CMZ/INSTALL=LOCATER,<files_directory>,[<option>,...]'
   mess ' '
   mess ' The files_directory must exist before calling cmz'
   mess ' '
   mess ' Compile and selection options '
   mess ' DEBUG    - debug option'
   mess ' LIB      - create source and library, delete object files (default)'
   mess ' SOURCE   - create source code only'
   mess ' OBJECT   - create source and object code only'
   mess ' ALL      - create source, object and library'
   mess ' INC      - Use cpp style #include commands for +SEQ'
   mess ' GNU      - Use the g77 and gcc compilers '
   mess ' CAR      - create car file and exit '
   mess ' CALIB    - Calibration code has been selected. '
   mess ' DBGPRNT  - Debug print lines activated '
   mess ' DBGPRNTX - Extreme Debug print lines activated '
   mess ' BXIMPROVE- BXIMPROVE selected '
   mess ' BXUSER   - BXUSER selected '
   mess ' NOCCODE  - Deselecting C code '
   mess ' HISTOG   - Fast Fuzzy Histograms Selected '
   mess ' LOCVERSION:x.xx/xx - select locater version x.xx/xx'
   mess ' CBOFFVERSION:x.xx/xx - select cboff version x.xx/xx'
   mess ' '
   mess 'example: cmz -install locater srccboff    (for Unix machines)'
   mess '         CMZ /INSTALL=LOCATER,[.SRCCBOFF] (for VMS machines)'
   mess '  '
   mess ' Warning: SOURCE and LIB may not be compatible on all machines.'
   mess ' '

 RETURN
+ADD,TCCOMMON,TC_MACRO,0.
*CMZU:  2.01/05 02/04/96  16.18.40  by  F.-H.Heinsius
*CMZU:  2.01/02 14/03/96  18.05.20  by  Rafik Ouared
+ADD,TCCOMMON,TC_MACRO,111.
*             14/03/96 R.Ouared
*             TCLIFT modified to include 2 arrays: MTSV1,MTSV2
+ADD,TCCOMMON,TC_MACRO,191.
      INTEGER         MTSV1(5),MTSV2(5)
+ADD,TCCOMMON,TC_MACRO,198.
     &               ,MTSV1,MTSV2
+REP,TCCOMMON,TC_MACRO,889-907.
*     some new parameters
* V3CHCT - TCVER3: 2-prong charge conservation flag:0=no cons.1=do cons.
* VHBD   - outlier rejection flag: 1-reject more outliers 0-reject normal outl
* VHOR   - outlier chi^2 threshold
* VHSP   - hit phi scaling factor
* VHSZ   - hit zed scaling factor
*
* The parameters that end in D are the Defaults
*
      INTEGER V3MNHD
      INTEGER V3CHCD,VHBDD
      REAL V3COND, V3DISD, V3BACD
      REAL VHMCHD,VHGCHD,VHCFLD,VHERCD,VHCTD
      REAL VHORD, VHSPD, VHSZD
      PARAMETER ( V3COND = 0.0, V3DISD=20.0, V3BACD=100.0 )
      PARAMETER ( V3MNHD = 7 )
      PARAMETER ( VHMCHD = 50.0, VHGCHD = 1.8, VHCFLD = 0.0 )
      PARAMETER ( VHERCD = 1.0, VHCTD = 0.001 )
      PARAMETER ( V3CHCD = 1, VHBDD = 1)
      PARAMETER ( VHORD =  25.0, VHSPD = 0.917, VHSZD = 0.886 )
*
      INTEGER VERTVR,VFITVR,V3MNHT
      REAL V3CONF, V3DIST, V3BACK
      REAL VHMCH,VHGCH,VHCFL,VHERC, VHCT
      INTEGER V3CHCT,VHBDT
      REAL VHORT,VHSPT,VHSZT
*
      COMMON /TCVTCM/
     &     VERTVR, VFITVR,
     &     V3CONF, V3DIST, V3BACK, V3MNHT,
     &     VHMCH, VHGCH, VHCFL, VHERC, VHCT,
     &     V3CHCT, VHBDT, VHORT, VHSPT, VHSZT

+ADD,TC_MAIN,CJINIT,0.
*CMZU:  2.01/10 22/01/97  15.09.56  by  M.T. Lakata
*CMZU:  2.01/05 02/04/96  14.52.15  by  F.-H.Heinsius
+ADD,TC_MAIN,CJINIT,358.
*     some more parameters from tcvtc2 common
      V3CHCT = V3CHCD
      VHBDT = VHBDD
      VHORT = VHORD
      VHSPT = VHSPD
      VHSZT = VHSZD
+ADD,TC_MAIN,CJINIT,529.
      CALL FFKEY('V3CH',V3CHCT, 1,'I') ! charge conservation flag
+ADD,TC_MAIN,CJINIT,535.
      CALL FFKEY('VHBD',VHBDT, 1,'I') ! outlier flag
      CALL FFKEY('VHOR',VHORT, 1,'R') ! outlier chi^2 cutoff
      CALL FFKEY('VHSP',VHSPT, 1,'R') ! hit error in rho, scaling factor
      CALL FFKEY('VHSZ',VHSZT, 1,'R') ! hit error in zed, scaling factor
+REP,TC_MAIN,CJINIT,549.
294     FORMAT(' Locater Card VERT ',I1,' is out range [1,3]')
+REP,TC_MAIN,CJINIT,555.
295     FORMAT(' Locater Card VFIT ',I1,' is out range [1,2]')
+ADD,TC_MAIN,TCINIT,0.
*CMZU:  2.01/13 14/03/97  15.46.32  by  M.T. Lakata
*CMZU:  2.01/12 25/02/97  10.45.51  by  Michael Doser
*CMZU:  2.01/09 15/11/96  17.39.10  by  Michael Doser
+ADD,TC_MAIN,TCINIT,167.
*                        15 Nov 1996, M. Doser
*                        Include the SVX calibration constants, both
*                        for global as for detector specific corrections
*
*                        10 December 1996 C.A.Meyer
*                        Forgot to declare ISTR and IDET
*
*                        12 February 1997 M.Doser
*                        move the calculation of the SVX geometry with the
*                        correct calibration constants from here to TCINIB
*                        (new routine), which is called from CBSRUN
*
*                        14 March 1997 M.Lakata
*                        removed the redundant code that is now in TCINIB

+REP,TC_MAIN,TCINIT,180-182.
      INTEGER        ISTR,IDET

*
      REAL           ANGLE,RADAN,RADIN,PHI,ALPRG,ALPLF
+REP,TC_MAIN,TCINIT,217.
      REAL    ASVXT(3)
+ADD,TC_MAIN,TCINIT,262.
&CDE,SVCALI.
+DEL,TC_MAIN,TCINIT,292.
+DEL,TC_MAIN,TCINIT,392-396.
+REP,TC_MAIN,TCINIT,625-629.
&SELF,IF=TVXT.
      WRITE(LLOG,*) 'TCINIT: ASVXTP(1,3)..',(ASVXTP(I),I=1,3)
+REP,TC_MAIN,TCINIT,787-849.
* default values for SVX, most likely changed with proper database values
*
      SVXANG = 0.366
      SVXGRD = 1.25
      SVXGDX = 0.
      SVXGDY = 0.
      SVXGDZ = 0.
      DO I=1,15
         SVXDDX(I) = 0.
         SVXDDY(I) = 0.
         SVXDDZ(I) = 0.
         SVXDDA(I) = 0.
         SVXDDB(I) = 0.
         SVXDDG(I) = 0.
      ENDDO

* positions are now calculated in TCINIB, which
* is also called by CBSRUN for new runs

      CALL TCINIB
+ADD,TC_MAIN,TCTRAK,0.
*CMZU:  2.01/14 07/08/97  00.10.38  by  Rafik Ouared
+ADD,TC_MAIN,TCTRAK,157.
*                       04/08/97  R.Ouared
*                       Add M.Doser's code to solve the following problem:
*                      the TCTR banks are still filled from the
*                      previous event ! garbage disposal necessary !
*                      Otherwise, this leads to a fatal error in GTSOFF,
*                      which uses the TTKS bank to look into the LTCTK
*                      bank
+REP,TC_MAIN,TCTRAK,659-660.
*        M.D. 6.8.97 : but the TCTR ... banks are still filled from the
*                      previous event ! garbage disposal necessary !
*                      Otherwise, this leads to a fatal error in GTSOFF,
*                      which uses the TTKS bank to look into the LTCTK
*                      bank
*
      IF(LGABT .OR. (IERR .NE. 0)) THEN
         IF(LTCTR.GT.0) THEN
            CALL MZDROP(IXSTOR,LTCTR,'L')
            LTCTR = 0
         ENDIF
         IF(LTCHX.GT.0) THEN
            CALL MZDROP(IXSTOR,LTCHX,'L')
            LTCHX = 0
         ENDIF
         IF(LTCVX.GT.0) THEN
            CALL MZDROP(IXSTOR,LTCVX,'L')
            LTCVX = 0
         ENDIF
         RETURN
      ENDIF
+ADD,TC_MAIN,ZBFORM,0.
*CMZU:  2.01/02 14/03/96  12.21.27  by  F.-H.Heinsius
*--  Obsolete
*CMZU:          14/03/96  18.07.23  by  Rafik Ouared
+DEL,TC_MAIN,ZBFORM,16-150.
+ADD,TC_MAIN,TCFORM,0.
*CMZU:  2.01/12 03/03/97  16.09.15  by  Rafik Ouared
*CMZU:  2.01/02 14/03/96  12.19.43  by  F.-H.Heinsius
+ADD,TC_MAIN,TCFORM,25.
*        Modified:       14,March,1996 R.Ouared
*                        NTSV1,NTSV2 arrays added for svtx in tsvpos.
*                        Bug in NRSVX data fixed up. NTVXT array
*                        updated to include 4th bank in TSVPOS after
*                        M.Doser request.
+ADD,TC_MAIN,TCFORM,46.
      INTEGER  NTSV1(5),NTSV2(5)
+REP,TC_MAIN,TCFORM,72-73.
      DATA    NTVXT  /4HTVXT, 4, 4,    6,2/
      DATA    NRSVX  /4HRSVX, 0, 0,    0,1/
      DATA    NTSV1  /4HTSV1, 0, 0,    0,0/
      DATA    NTSV2  /4HTSV2, 0, 0,    0,0/
+ADD,TC_MAIN,TCFORM,117.
         MTSV1(I) = NTSV1(I)
         MTSV2(I) = NTSV2(I)
+REP,TC_MAIN,TCFORM,147.
*      CALL MZFORM('RSVX','/2I 8F',     MRSVX(5))
      CALL MZFORM('RSVX','-B',     MRSVX(5))
      CALL MZFORM('TSV1','/2I',        MTSV1(5))
      CALL MZFORM('TSV2','/3I 12F 1I', MTSV2(5))
+ADD,TC_MAIN,TCINIB,*.
&DECK,TCINIB.
*CMZU:  2.01/12 12/02/97  16.04.52  by  Michael Doser
*-- Author : Michael Doser
      SUBROUTINE TCINIB
*
*        Author:         Michael Doser
*
*        Creation date:  12 February  1997.
*
*        Description:    This routine is called once at the beginning
*                        of each run. It's purpose is to calculate the
*                        positions of the 1920 strips of the SVX
*
*        Subroutines Referenced:
*
*        Edit History:
*
*     **************************************************************************
*
&SEQ,IMPNONE.
*
*        Internal variables:
*
      INTEGER        I,II
      INTEGER        ISTR,IDET
 
*
      REAL           ANGLE,RADAN,RADIN,PHI,ALPRG,ALPLF
      REAL           SECANG,ANGLS,PHIA,COSANG,SINANG
*
&CDE,CBUNIT.
*
&CDE,PI2PI.
*
&CDE,TPPRMS.
&CDE,SVCALI.
*
*
*        ########################################
*        # geometry for silicon vertex detector #
*        ########################################
*
      SECANG= PI2/15.000
      RADIN = PI/180.000
*
      WRITE (LLOG,*) ' Updating SV geometry with database corrections'
      DO 8000 I = 0,1919
*
* find x,y and radius for strips of first module; start with
* start with strip 128, calculate the others from there
*
        ISTR=MOD(I,128)
        IDET=I/128+1
        II=128*(IDET-1)+(127-ISTR)
        IF(ISTR.EQ.0) THEN
           XSVXTP(II)=SVXGRD*COS(SVXANG*RADIN-SECANG*(IDET-1)) + SVXGDX
           YSVXTP(II)=SVXGRD*SIN(SVXANG*RADIN-SECANG*(IDET-1)) + SVXGDY
           IF(XSVXTP(II).EQ.0.) XSVXTP(II)=1.0E-06
           PHIA = ATAN2(YSVXTP(II),XSVXTP(II))
           IF(PHIA.LT.0.) PHIA=PHIA+PI2
           ANGLS = PHIA+PIHLF-ASVXTP(3)+SVXDDG(IDET)
           COSANG=COS(ANGLS)
           SINANG=SIN(ANGLS)
 
           XSVXTP(II) = XSVXTP(II) + SVXDDX(IDET)
           YSVXTP(II) = YSVXTP(II) + SVXDDY(IDET)
*
*  strips 1..127
        ELSE
*
           XSVXTP(II)=XSVXTP(128*(IDET-1)+127)
     &                  + ASVXTP(1)*FLOAT(ISTR)*COSANG
           YSVXTP(II)=YSVXTP(128*(IDET-1)+127)
     &                  + ASVXTP(1)*FLOAT(ISTR)*SINANG
 
        ENDIF
 
        PHI = ATAN2(YSVXTP(II),XSVXTP(II))
        IF(PHI.LT.0.) PHI=PHI+PI2
        IF (II.LT.128) THEN
           RSVXTP(II)=SQRT(XSVXTP(II)*XSVXTP(II)+
     &                       YSVXTP(II)*YSVXTP(II))
        ENDIF
        FSVXTP(II)=PHI
 
&SELF,IF=TCINIT.
        IF(I.LT.128) THEN
          WRITE(LLOG,*) 'TCINIT: STRIP ..I,II..',I,II
          WRITE(LLOG,*) 'TCINIT: XSVXTP,YSVXTP,FSVXTP..'
          WRITE(LLOG,*)  XSVXTP(II),YSVXTP(II),FSVXTP(II)
        ELSE
          WRITE(LLOG,*) 'TCINIT: STRIP ..I..',I
          WRITE(LLOG,*) 'TCINIT: XSVXTP,YSVXTP,FSVXTP..'
          WRITE(LLOG,*)  XSVXTP(I),YSVXTP(I),FSVXTP(I)
        ENDIF
&SELF.
*
 8000  CONTINUE
*
      RETURN
      END
 
+ADD,TC_RAWS,TJDCGT,0.
*CMZU:  2.01/08 21/06/96  14.24.50  by  F.-H.Heinsius
+ADD,TC_RAWS,TJDCGT,501.
*
*        For Monte Carlo Data add the global time_zero jitter to the hit.
*
                  IF(LTYPMC.AND.SMEATC)
     $                 Q(ITJDC+5) = MAX(0.0001,Q(ITJDC+5)+TIM0MC)

+REP,TC_RAWS,TJDCGT,513-518.
*                    Changed 19.6.96: smear distances, no more call
*                    to TJTIME. See comments in TSMUSH.
*
                  IF(LTYPMC.AND.SMEATC) THEN
                     CALL TSMUSH(LYR,ISEC,ITJDC)
CFHH                 ISECTJ = ISEC
CFHH                 ITJDTJ = ITJDC
CFHH                 CALL TJTIME
+ADD,TC_RAWS,TJDCGT,642.
*
*        For Monte Carlo Data add the global time_zero jitter to the hit.
*
                  IF(LTYPMC.AND.SMEATC)
     $                 Q(ITJDC+5) = MAX(0.0001,Q(ITJDC+5)+TIM0MC)

+REP,TC_RAWS,TJDCGT,651-656.
*                    Changed 19.6.96: smear distances, no more call
*                    to TJTIME. See comments in TSMUSH.
*
                     IF(LTYPMC.AND.SMEATC) THEN
                        CALL TSMUSH(LYR,ISEC,ITJDC)
CFHH                    ISECTJ = ISEC
CFHH                    ITJDTJ = ITJDC
CFHH                    CALL TJTIME
+ADD,TC_RAWS,TJTIME,0.
*CMZU:  2.01/08 26/06/96  16.59.32  by  F.-H.Heinsius
+REP,TC_RAWS,TJTIME,722-723.
CFHH* FHH 6/96: not for MC
      IF (.NOT. OLDJDC .AND. (.NOT.LTYPMC)) THEN
+REP,TC_RAWS,TJTIME,819.
* FHH 6/96: not for MC
      IF (.NOT.OLDJDC .AND. (.NOT.LTYPMC)) THEN
+REP,TC_RAWS,TJTIME,925.
* FHH 6/96: not for MC
      IF (OLDJDC.AND.(ISECTJ.EQ.1).AND.(.NOT.LTYPMC)) THEN
+ADD,TC_RAWS,TPWPOS,0.
*CMZU:  2.01/05 02/04/96  15.58.55  by  F.-H.Heinsius
*CMZ :  2.01/04 18/03/96  15.37.30  by  F.-H.Heinsius
+ADD,TC_RAWS,TPWPOS,53.
*                       18 March 1996---F.H.Heinsius
*                       Error DPHI was negative because wire distance
*                       APWCTP is negative!
*
+REP,TC_RAWS,TPWPOS,238.
            DPHI = 1.50*ABS(APWCTP(1))
+REP,TC_RAWS,TPWPOS,257.
            DPHI= 1.50*ABS(APWCTP(2))
+REP,TC_RAWS,TPWPOS,265-266.
         DX   = ABS(YPWC)*DPHI
         DY   = ABS(XPWC)*DPHI
+REP,TC_RAWS,TPWPOS,275.
*MTL.  Emperical scaling factor = 3.1. Perhaps might be 3.0
*MTL.  but 3.1 works well. (gives pull width of 1.009)
         Q(LTPCH+10) = DPHI/3.1
+ADD,TC_RAWS,TSMINI,0.
*CMZU:  2.01/08 26/06/96  17.00.21  by  F.-H.Heinsius
+ADD,TC_RAWS,TSMINI,24.
*        24.6.96: F.-H.Heinsius
*           Change global time offset from -5.5ns to +2.5ns
*
+REP,TC_RAWS,TSMINI,168.
CFHH     TIM0MC = -0.0040 - 0.0030*TIM0MC
         TIM0MC =  0.0040 - 0.0030*TIM0MC
+ADD,TC_RAWS,TSMUSH,0.
*CMZU:  2.01/08 19/06/96  16.00.04  by  F.-H.Heinsius
+ADD,TC_RAWS,TSMUSH,35.
*        19.6.96: F.-H.Heinsius: Instead of calculating a drifttime from
*                 the smeared positions modify positions directly. Do
*                 not call TJTIME again to recalculate positions.
*                 The corrections were originally calculated based on
*                 position not on time. However, this leads to an
*                 inconsistency because the TIME value does not reflect
*                 the corresponding position. This should not be a
*                 problem because it used only in the obsolete TCPATT
*                 subroutines. But nevertheless be careful!
*
*                 TIME is not updated correctly, use ONLY POSITIONS after this
*                 Subroutine!
*
*                 Also fix bug in relative distance to wire for new JDC.
+ADD,TC_RAWS,TSMUSH,76.
&SEQ,XYSENS.
&SEQ,TCFLAG.
+REP,TC_RAWS,TSMUSH,92-95.
*        Get the wire position (fixed for new JDC, 19.6.96)
*
      IF( (.NOT.OLDJDC).AND.(LAYER.LT.6) ) THEN
*
         XWIRS = XSENSE(LAYER)
         YWIRS = YSENSE(LAYER)
*
      ELSE
*
         XWIRS = RLYRTC(LAYER)
         YWIRS = Q(LTJST + (ISEC-1)*23 + LAYER)
*
      ENDIF
+REP,TC_RAWS,TSMUSH,102-106.
*        a distance smearing in centimeters.

      DIST = TSRPHI(LAYER, R )

      IF (R.GT. 1E-10) THEN
         DISTX = DIST * (Q(KTJDC+8)-XWIRS)/R
         DISTY = DIST * (Q(KTJDC+9)-YWIRS)/R
      ELSE
         DISTX = 0
         DISTY = DIST
      END IF

      Q(KTJDC+ 8) = Q(KTJDC+ 8)+DISTX
      Q(KTJDC+ 9) = Q(KTJDC+ 9)+DISTY
      Q(KTJDC+10) = Q(KTJDC+10)-DISTX
      Q(KTJDC+11) = Q(KTJDC+11)-DISTY
*
*        Then multiply this by (1/v_drift) in cm/micro-sec to update time.
*
+REP,TC_RAWS,TSMUSH,119-121.
*        Add the global time_zero jitter to the hit. (done in TJDCGT)
*
**      Q(KTJDC+5) = MAX(0.0001,Q(KTJDC+5)+TIM0MC)
*

+ADD,TC_RAWS,TSMUSH,149.

*     update Z Position:

      Q(KTJDC+12) = Q(KTJDC+12)+DELZ


+ADD,TC_RAWS,TSRPHI,0.
*CMZU:  2.01/05 02/04/96  19.16.58  by  F.-H.Heinsius
+REP,TC_RAWS,TSRPHI,39-40.
      REAL E,D,R,S,C(6),TL,AR
      INTEGER LAYER,B,I
+ADD,TC_FFUZ,TCCALB,0.
*CMZU:  2.01/14 03/06/97  12.03.23  by  Matthias Heinzelmann
+ADD,TC_FFUZ,TCCALB,19.
*		       03/06/97  Matthias Heinzelmann
*		       proper initialization of IPATT,RPATT,DPATT,
*		       LINK_HIT
*                      changed loopindex from MCLUSTER to MSECTOR
*		       -> removes dependence of the reconstruction
*		          from previous event
+ADD,TC_FFUZ,TCCALB,137.

* MH >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
         DO I=1,MHIT
           DO J=1,MIINF
             IPATT(I,J)=0
           ENDDO
         ENDDO
* MH <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<

+REP,TC_FFUZ,TCCALB,210-216.

* MH >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
            DO I=1,MHIT
              DO J=1,MRINF
                RPATT(I,J) = 0.
              ENDDO
            END DO
* MH <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<

         DO IHIT  = 1,NR_ALL_HITS
+ADD,TC_FFUZ,TCCALB,244.

* MH >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
         DO I=1,MHIT
           DO J=1,MDINF
             DPATT(I,J)=0.
           ENDDO
         ENDDO
* MH <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<

+REP,TC_FFUZ,TCCALB,287.
* MH >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
*        DO ISECTOR = 1,MCLUSTER
         DO ISECTOR = 1,MSECTOR
* MH <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
+ADD,TC_FFUZ,TCCALB,294.

* MH >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
         DO J=1,MHIT
           LINK_HIT(J)=0
         ENDDO
* MH <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<

+ADD,TC_FFUZ,TCFCIRD,0.
*CMZU:  2.01/14 03/06/97  11.40.32  by  Matthias Heinzelmann
+ADD,TC_FFUZ,TCFCIRD,26.
*                       03/06/97  Matthias Heinzelmann
*			initialize NHITNW, RFFIT
*			-> removes dependece of reconstruction
*			   from previous event
*
+REP,TC_FFUZ,TCFCIRD,45.
      INTEGER ISTEP,NSTEP,I,NTRY,J
+ADD,TC_FFUZ,TCFCIRD,115.

* MH >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
      NHITNW = 0
      DO I=1,MFIT
        DO J=1,MFIT_INF
          RFFIT(I,J)=0.
        ENDDO
      ENDDO
* MH <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<

+ADD,TC_FFUZ,TCFFUZ,0.
*CMZU:  2.01/14 03/06/97  11.38.30  by  Matthias Heinzelmann
+ADD,TC_FFUZ,TCFFUZ,38.
*         	       03/06/97 Matthias Heinzelmann
*                      call TCKON1(0) at the beginning of the
*                      loop (initialization)
*
+ADD,TC_FFUZ,TCFFUZ,179.

* MH >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
         CALL TCKON1(0)
* MH <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<

+REP,TC_FFUZ,TCFFUZ,259.
* MH >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
*           CALL TCKON1(0)
* MH <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
+ADD,TC_RADON1,RADON1,0.
*CMZU:  2.01/05 02/04/96  18.51.11  by  F.-H.Heinsius
+DEL,TC_RADON1,RADON1,52.
+ADD,TC_RADON1,TCFILL,0.
*CMZU:  2.01/02 14/03/96  15.31.57  by  F.-H.Heinsius
+REP,TC_RADON1,TCFILL,95-96.
         END DO
         NR_TRACKS_FOUND = ITRACK
+ADD,TC_CIRC,TCCIRC,0.
*CMZU:  2.01/09 01/10/96  16.02.54  by  Rafik Ouared
+ADD,TC_CIRC,TCCIRC,98.
*                 **>   01 october 1996 --- R.Ouared
*                       Max number of hit per track is now 50.
*
+ADD,TC_CIRC,TCCIRC,234.
***more than 50 hits per track not allowed.
*
         IF(N.GT.50) THEN
            N=MIN(N,50)
            WRITE(LLOG,*) 'TCCIRC: MORE THAN 50 HITS PER TRACK FOUND...'
         ENDIF
*
+ADD,TC_CIRC,TXTGET,0.
*CMZU:  2.01/07 17/06/96  17.35.28  by  F.-H.Heinsius
*CMZU:  2.01/05 02/04/96  16.10.05  by  F.-H.Heinsius
*CMZ :  2.01/01 13/03/96  15.14.15  by  F.-H.Heinsius
+ADD,TC_CIRC,TXTGET,46.
*                       13 Mar 1996, F.-H.Heinsius
*                       Return error instead of error**2 for PWC/SVX.
*
*                       17 March 1996 M.Lakata
*                       Some array values (SR,SP) are now filled for
*                       pwc and svx.
+REP,TC_CIRC,TXTGET,104-108.
* ITPCL : pointer to pwc bank
* J     : loop index
* ITVTX : pointer to svx bank

       INTEGER IPWC,IJDC,ITPCL,J,ITVTX
+REP,TC_CIRC,TXTGET,153.
            ITPCL = LQ(LTPWC-I)
+REP,TC_CIRC,TXTGET,159-173.
               IF(IQ(ITPCL+1).EQ.ITRK) THEN
*
*                    Add this hit to the track:
*
                  NPTS = NPTS + 1
                  X1(NPTS) = Q(ITPCL+5)
                  Y1(NPTS) = Q(ITPCL+6)
                  Z1(NPTS) = Q(ITPCL+7)
*
*                    Store the errors in (x,y,z). Note that
*                    the error is z is +/- 25cm 1 sigma!
*
                  D1(NPTS) = Q(ITPCL+10)
                  D2(NPTS) = Q(ITPCL+11)
                  D3(NPTS) = 25.0

*NEW
                  R1(NPTS)   = Q(ITPCL+8)
                  ILAY(NPTS) = 0
                  ISC(NPTS)  = 0
                  IRS(NPTS)  = 0
                  TIM(NPTS)  = 0.0
                  AMPL(NPTS) = 1.0
                  AMPR(NPTS) = 1.0
* guess at the pwc wire thickness, 8 microns??
                  SR(NPTS)   = 0.0008
* there is a bug in TPWPOS that makes the error negative:
                  SP(NPTS)   = ABS(Q(ITPCL+12))
*ENDNEW
+REP,TC_CIRC,TXTGET,182.
               ITPCL = ITPCL + LENCL
+REP,TC_CIRC,TXTGET,193.
         ITVTX = LQ(LTVXT-2)
+REP,TC_CIRC,TXTGET,200-216.
            IF(IQ(ITVTX+1).EQ.ITRK) THEN
*
*                    Add this hit to the track:
*
               NPTS = NPTS + 1
               X1(NPTS) = Q(ITVTX+4)
               Y1(NPTS) = Q(ITVTX+5)
               Z1(NPTS) = Q(ITVTX+6)
*
*                    Store the errors in (x,y,z). Note that
*                    the error is z is +/- 3.5cm 1 sigma!
*
               D1(NPTS) = Q(ITVTX+10)
               D2(NPTS) = Q(ITVTX+11)
               D3(NPTS) = 3.5
*NEW
*         PHI  =Q(ITVTX+8)
*         PHIP =Q(ITVTX+9)
*         DXX  =Q(ITVTX+10)
*         DYY  =Q(ITVTX+11)
*         DXY  =Q(ITVTX+12)
*         DPHI =Q(ITVTX+13)
*         DRAD =Q(ITVTX+14)
*         DEDX =Q(ITVTX+15)
*         BPAMP=IQ(ITVTX+16)

               R1(NPTS)   = Q(ITVTX+7)
* bogus values
               ILAY(NPTS) = 0
               ISC(NPTS)  = 0
               IRS(NPTS)  = 0
               TIM(NPTS)  = 0.0
               AMPL(NPTS) = 1.0
               AMPR(NPTS) = 1.0
*     guess at the svx thickness, 8 microns??
               SR(NPTS)   =Q(ITVTX+14)
               SP(NPTS)   =Q(ITVTX+13)
*ENDNEW
            ENDIF
*
            ITVTX = ITVTX + LENCLS
+DEL,TC_CIRC,TXTGET,268.
+ADD,TC_HELX,TCHXLD,0.
*CMZU:  2.01/03 15/03/96  12.51.46  by  F.-H.Heinsius
+REP,TC_HELX,TCHXLD,478-482.
*                    the error is z is +/- 3.5cm 1 sigma!
*
               Q(ITCHX+4) = Q(ITVXT+10)**2
               Q(ITCHX+5) = Q(ITVXT+11)**2
               Q(ITCHX+6) = 3.5**2
+ADD,TC_HELX,TPCNCT,0.
*CMZU:  2.01/05 02/04/96  16.21.53  by  F.-H.Heinsius
+REP,TC_HELX,TPCNCT,274.
      INTEGER    ITRK,ITCTK,ITPWC,KTPCL,NRHITS,NRJDC
+REP,TC_HELX,TPCNCT,554.
*        create a cluster bank for this pwc. KTPCL will point into
+REP,TC_HELX,TPCNCT,561.
         CALL MZLIFT(IXTPWC,KTPCL,LTPWC,ILINK,MTPCL,0)
+REP,TC_HELX,TPCNCT,573.
            IQ(KTPCL+1) = -1
+REP,TC_HELX,TPCNCT,582-583.
            IQ(KTPCL+2) = IQ(ITPWC+2)
            IQ(KTPCL+3) = CLRW(J,M)
+REP,TC_HELX,TPCNCT,590.
            Q(KTPCL+4)  = FLOAT(MAXWYR(M)-MWOF1) - CLNR(J,M)
+REP,TC_HELX,TPCNCT,602-610.
            Q(KTPCL+8) = RPWCTP(M)
            Q(KTPCL+9) = CLANG(J,M)
*
*              Now compute the (x,y) of the cluster as:
*              x = r * cos(phi)   y = r * sin(phi)
*
            Q(KTPCL+5) = RPWCTP(M) * COS(CLANG(J,M))
            Q(KTPCL+6) = RPWCTP(M) * SIN(CLANG(J,M))
            Q(KTPCL+7) = 0.0
+REP,TC_HELX,TPCNCT,617-630.
            Q(KTPCL+12) = Q(ITPWC+10)/SQCLR
*
*              The error in x and y can be computed simply as:
*              dx = y * dphi     dy = x * dphi
*
            Q(KTPCL+10) = Q(KTPCL+6) * Q(KTPCL+12)
            Q(KTPCL+11) = Q(KTPCL+5) * Q(KTPCL+12)
*
*              Advance the KTPCL pointer to allow us to store the next
*              hit.
*
&SEQ,TPCNCT23,T=PASS,IF=DEBHELIX.
*
            KTPCL = KTPCL + LENCL
+REP,TC_HELX,TPCNCT,979-983.
            KTPCL = LQ(LTPWC-M-2) + LENCL * (ICL-1)
*
*              Now stamp the cluster with the JDC track number.
*
            IQ(KTPCL+1) = I
+ADD,TC_VERT,TCVRTX,0.
*CMZU:  2.01/13 14/03/97  16.24.30  by  M.T. Lakata
+ADD,TC_VERT,TCVRTX,158.
*                       14 March 1997, M. Lakata
*                       fix some mzlift garbage collection errors
*
+REP,TC_VERT,TCVRTX,184.
*              or (if negative) is the index in XXTCVT
+ADD,TC_VERT,TCVRTX,311.
&SEQ,TCV3CM.
+REP,TC_VERT,TCVRTX,323.
        ITCVT = XXTCVT(-IVRTX)
+REP,TC_VERT,TCVRTX,974-975.
      IF (IVRTX.LT.0) THEN
         ITCVT = XXTCVT(-IVRTX)
      ELSE
         ITCVT = LQ(LTCVX-IVRTX)
      ENDIF

      CALL MZLIFT(IXTCVX,ITCVP,ITCVT,-1,MTCVP,0)
      IF (IVRTX.LT.0) THEN
         ITCVT = XXTCVT(-IVRTX)
      ELSE
         ITCVT = LQ(LTCVX-IVRTX)
      ENDIF
      ITCVP = LQ(ITCVT-1)
+REP,TC_VERT,TCVRTX,983-984.
* extra steps necessary in case of garbage collection
      IF (IVRTX.LT.0) THEN
         ITCVT = XXTCVT(-IVRTX)
      ELSE
         ITCVT = LQ(LTCVX-IVRTX)
      ENDIF
      ITCVP = LQ(ITCVT-1)
      ITCVH = LQ(ITCVP)
+ADD,TC_VERT,TCVER3,0.
*CMZU:  2.01/13 14/03/97  18.12.00  by  M.T. Lakata
*CMZU:  2.01/10 22/01/97  17.32.03  by  M.T. Lakata
*CMZU:  2.01/07 17/06/96  18.05.36  by  F.-H.Heinsius
*CMZU:  2.01/05 02/04/96  15.53.39  by  F.-H.Heinsius
+REP,TC_VERT,TCVER3,12-20.
*     rewritten      27 March 1996
*     reference      CBNote 285
*
*     March 14, 1997 XXTCVT is now a common zebra link bank, so that
*     garbage collection doesn't ruin the temporary banks (a problem
*     with large ntuples)
*
*     Goal: Event with lowest vertex multiplicity, each with
*     "good" chi^2.
*
*     Algorithm:
*     1. make "dummy" 1-prong vertices
*     2. combine pairwise all current vertices into composite vertices
*        if possible.
*        If no pairs possible, quit
*     3. construct all configurations out of composite vertices and
*        unused old vertices
*     4. Find configuration with lowest chi^2
*     5. replace current list of vertices with best configuration
*     6. goto 2
*
*     Note that this algorithm is much more robust than the previous one
*     from August 1995.  The algorithm has a well defined method and goal,
*     that is vertex multiplicity minimization.
*
*     There are a few speed ups.
*     1. in the first stage, only +/- pairs are combined.
*     2. there are some cuts to decide whether to merge vertices
*        based on some geometry cuts.  Currently, it is suggested not
*        to depend on these values.
*        it is suggested that v3di and v3cl be very liberal,
*        for example: v3di = 6.0cm and v3cl = 0.0
+REP,TC_VERT,TCVER3,30-165.
*
&SEQ,IMPNONE.
&SEQ,CBLINK.
&SEQ,TCLIFT.
&SEQ,TCSTAT.
&SEQ,CBUNIT.
&SEQ,TRKPRM.
&SEQ,TCVTCM.
&SEQ,CBSTOP.
&SEQ,TCV3CM.
*
*     CONSTANTS
*
*     MAXCNF     maximum number of configurations
*     MAXTRK     maximum number of tracks

      INTEGER MAXTRK,MAXCNF
      PARAMETER (MAXTRK=8,MAXCNF=80)

*     Internal variables:
*
*     INTRODUCTION ----------------------------------------------
*     FIRST      true for first time
      LOGICAL FIRST
      SAVE FIRST

*     IVRTX      number of all vertices
*     NTRAK      number of helix tracks

      INTEGER   NTRAK,IVRTX

*     PART -1 ---------------------------------------------------
*     VALIDT     array of tracks that meet minimum requirements
      LOGICAL VALIDT(MAXTRK)

*     PART I   --------------------------------------------------
*     FULMSK     mask of all tracks "on"
*     NTRAKS     number of good tracks
*     ITRAK      track index
*     ITCTR      is the pointer to the TCTR bank
*     NHITS      number of hits on a track
*     IERR       error of vertex fit
*     CIRERR     circle error
*     HLXERR     helix error
*     J,K        loop variables
*     NTR        number of tracks to try for fit
*     TRKS       list of tracks to try
*     ITCHX1     temp pointer to tchx
*     NDROP      number of tracks dropped by tcvrtx
*     IDROP      dummies array for tcvrtx
*     TRKSVT     number of trks for a vertex,
*     VON        flag: 0=vertex off, 1=vertex used,2=vertex proposed
*     CHARGE     charge of vertex
*     MASK       bit-mask for each vertex. bit N is set for track N
*     OLDVRT     current number of "old" vertices
*     XGUES  is a guess to the vertex position.

      INTEGER FULMSK,NTRAKS,ITRAK,ITCTR,NHITS,IERR
      INTEGER CIRERR,HLXERR,J,K,NTR,TRKS(MAXTRK)
      INTEGER ITCHX1
      INTEGER NDROP,IDROP(MAXTRK)
      INTEGER TRKSVT(MAXVRT),VON(MAXVRT),CHARGE(MAXVRT)
      INTEGER MASK(MAXVRT),OLDVRT
      REAL    XGUES(3)

*     PART II  --------------------------------------------------
*     NEWVRT     current number of "new" (composite) vertices
*     IV,IV1,IV2 vertex indices
*     IVRTX1     similar to ivrtx
*     IVRTX2     similar to ivrtx
*     ITCVT      pointer to a TCVT bank (vertex parameters)
*     IFAIL      flag from RINV
*     TEMP       temp space
*     DX,DY,DZ   distance from vertex to vertex
*     R          distance between two vertices
*     P          probability/confidence level
*     ERR        3x3 error matrix
*     CHI        chi2 for a vertex (from TCVRTX)

      INTEGER NEWVRT
      INTEGER IV,IV1,IV2,IVRTX1,IVRTX2
      INTEGER ITCVT1,ITCVT2
      INTEGER IFAIL,TEMP(3)
      REAL    DX,DY,DZ
      REAL    R,P
      REAL    ERR(3,3),CHI

*     PROB    returns confidence level from chi2 and ndf
      REAL     PROB
      EXTERNAL PROB

*     PART III --------------------------------------------------
*     MAXNV   maximum number of desired vertices
*     NCONF   number of found valid configurations
*     NVLIST  temporary vertex list for a configuration
*     CUMASK  temporary bit-mask of all tracks
*     NVCONF  number of vertices in a configuration
*     COLIST  vertex list for a configuration
*     VLIST   temporary list of vertices in a configuration
*     NV      number of vertices in each configuration

      INTEGER NVLIST
      INTEGER NCONF,CUMASK,NVCONF(MAXCNF),NV
      INTEGER COLIST(MAXTRK,MAXCNF),VLIST(MAXTRK),MAXNV

*     PART IV  --------------------------------------------------
*     BESTCF  index of best configuration
*     NDF     number of degrees of freedom
*     ITCVT   temp pointer to tcvt
*     ICONF      configuration index
*     ITCVP  is a pointer to a TCVP bank (track parameters)
*     BESTCH  best Chi2 for a configurations
*     X,Y,Z   primary vertex position
*     PX,PY,PZ vertex momentum

      INTEGER BESTCF,ITCVT,NDF
      INTEGER ICONF, ITCVP
      REAL    BESTCH
      REAL    X,Y,Z,PX,PY,PZ

*     PART V   --------------------------------------------------
*     IGOOD   final count of good vertices
      INTEGER IGOOD

*     PART VI  --------------------------------------------------
*     ITCVH  is a pointer to a TCVH bank (helix parameters)
      INTEGER ITCVH
*
*     FOR ERRORS
+REP,TC_VERT,TCVER3,171-180.
      DATA IERRL /10*0/
      DATA FIRST /.TRUE./
      DATA NDROP /0/
*
***************************************************************************
*     INTRODUCTION
***************************************************************************

+REP,TC_VERT,TCVER3,190-213.
*     If any changes to defaults, print out warning, otherwise
*     WELCOME to TCVER3 !

         WRITE(LLOG,294)
 294     FORMAT(' *** TCVER3 Initialized')
         CALL TCVPCR('V3CL',V3CONF,V3COND)
         CALL TCVPCR('V3DI',V3DIST,V3DISD)
         CALL TCVPCR('V3BA',V3BACK,V3BACD)
         CALL TCVPCI('V3MH',V3MNHT,V3MNHD)
         CALL TCVPCI('V3CH',V3CHCT,V3CHCD)
         CALL MZLINK(IXTCVX,'/TCV3CM/',
     &     XXTCVT(1),XXTCVT(1),XXTCVT(MAXVRT))
+REP,TC_VERT,TCVER3,229-241.



***************************************************************************
* PART -1
* sort out good tracks/bad tracks
***************************************************************************
      NTRAKS = 0

*
*     First, we assign "1-prong" vertices to each "good" track
*     and make a list of vertices from thse
*
      DO 10 ITRAK = 1,NTRAK
*
*     Get the address for this track's TCTR bank.
*
         ITCTR = LQ(LTCTR-ITRAK)
*
*     Make sure that this track has not been connected to a vertex.
+REP,TC_VERT,TCVER3,267-350.
            NTRAKS = NTRAKS + 1
            VALIDT(ITRAK) = .TRUE.
         ELSE
            VALIDT(ITRAK) = .FALSE.
         ENDIF
 10   CONTINUE

***************************************************************************
* part 0
* try fitting all tracks together. If this fails, then try another method
* Since most events have no secondary vertices, this should save some CPU
* time.
***************************************************************************

      NTR = 0
      DO 20 ITRAK=1,NTRAK
         IF(VALIDT(ITRAK)) THEN
            NTR = NTR+1
            TRKS(NTR) = ITRAK
         ENDIF
 20   CONTINUE

*     Create a vertex and momentum bank for this vertex.
      IVRTX = 1
      CALL MZLIFT(IXTCVX,XXTCVT(IVRTX),XXTCVT(IVRTX),1,MTCVT,0)
*     call desired fitter
      IF(VFITVR.EQ.1) THEN
         CALL TCVRTX(NTR,TRKS,-IVRTX,NDROP,IDROP,
     &        IERR,XGUES)
      ELSEIF (VFITVR.EQ.2) THEN
         CALL TCVRHX(NTR,TRKS,-IVRTX,IERR,XGUES)
      ENDIF

************* CUT3 ************
      IF(IERR.GT.100.OR.NDROP.GT.0) THEN
*     If this was bad, then drop the fit, and forget it ever happened
*     continue on with bottom-up approach

         CALL MZDROP(IXSTOR,XXTCVT(IVRTX),'L')
         IVRTX=IVRTX-1
      ELSE
*     hoorah, we got a good one on the first try
*     turn it on
         VON(IVRTX)    = 1
*     save it and quit!
         GOTO 4000
      ENDIF

***************************************************************************
*     PART 1
***************************************************************************

      FULMSK = 0
      NTRAKS = 0


*
*     First, we assign "1-prong" vertices to each "good" track
*     and make a list of vertices from thse
*
      DO 100 ITRAK = 1,NTRAK
*
*     Get the address for this track's TCTR bank.
*
         ITCTR = LQ(LTCTR-ITRAK)
*
*     Make sure that this track has not been connected to a vertex.
*
         IQ(ITCTR+3) = 0

         IF (VALIDT(ITRAK)) THEN
            TRKS(1) = ITRAK

* initial guess is average of first hit
            ITCHX1 = LQ(LTCHX-TRKS(1))
            DO 755 K=1,3
               XGUES(K)=Q(ITCHX1+K)
 755        CONTINUE

*     Create a vertex and momentum bank for this vertex.

            IVRTX = IVRTX + 1
            CALL MZLIFT(IXTCVX,XXTCVT(IVRTX),XXTCVT(IVRTX),1,MTCVT,0)
*     call desired fitter
            IF(VFITVR.EQ.1) THEN
               CALL TCVRTX(1,TRKS,-IVRTX,NDROP,IDROP,
     &              IERR,XGUES)
            ELSEIF (VFITVR.EQ.2) THEN
               CALL TCVRHX(1,TRKS,-IVRTX,IERR,XGUES)
            ENDIF
            IF(IERR.NE.100) THEN
               CALL MZDROP(IXSTOR,XXTCVT(IVRTX),'L')
               IVRTX=IVRTX-1
            ELSE
*     keep the vertex
               CHARGE(IVRTX) = NINT(Q(ITCTR+7))
*     number of tracks at this vertex
               TRKSVT(IVRTX) = 1
               VON(IVRTX) = 1
*     construct a mask (i know that we waste a bit, but that's fortran)
*     (itrak = 1...N)
               MASK(IVRTX) = ISHFT(1,ITRAK)
*     WE make a "FullMask" that contains all the tracks we
*     want.  When the configuration == the fullmask, then
*     all the tracks have been used.
               FULMSK = FULMSK+MASK(IVRTX)
               NTRAKS = NTRAKS + 1
            ENDIF
         ENDIF
 100  CONTINUE
      OLDVRT = NTRAKS
*
***********************************************************
* PART II
***********************************************************
*
*     new iteration
 1000 CONTINUE
*
*
      NEWVRT = 0
*     save the original number of vertices
      IVRTX1 = IVRTX-1
      IVRTX2 = IVRTX
*     and use it
      DO 200 IV1=1,IVRTX1
         IF (VON(IV1).EQ.0) GOTO 200
         DO 250 IV2=IV1+1,IVRTX2
            IF (VON(IV2).EQ.0) GOTO 250
*     try to merge these two, IV1 and IV2
*     if they are single prong and not neutral, don't merge
            IF (V3CHCT.GT.0.AND.TRKSVT(IV1).EQ.1.AND.TRKSVT(IV2).EQ.1)
     &           THEN
               IF (CHARGE(IV1)+CHARGE(IV2).NE.0) GOTO 250
            ENDIF


*     make some fast cuts first:
            ITCVT1 = XXTCVT(IV1)
            ITCVT2 = XXTCVT(IV2)
            DX = Q(ITCVT1+5)-Q(ITCVT2+5)
            DY = Q(ITCVT1+6)-Q(ITCVT2+6)
            DZ = Q(ITCVT1+7)-Q(ITCVT2+7)
            R = SQRT(DX*DX+DY*DY+DZ*DZ)

************ CUT1 **************
*     first check to see that distance between vertices is even close
C MTL changed on 22/1/97: Don't check distance for two 1prong vertices
C This allows one to still use VFIT 1 in the steering cards.
             IF(R.GT.V3DIST.AND..NOT.
     &           (TRKSVT(IV1).EQ.1.AND.TRKSVT(IV2).EQ.1)) GOTO 250

************ CUT2 **************
*     check confidence level for hypothesis that the two
*     vertices are the same
*     DEBUG: removed the offdiagonal elements....
            IF(V3CONF.GT.0.0) THEN
               ERR(1,1) = Q(ITCVT1+8 ) + Q(ITCVT2+8 )
               ERR(2,1) = Q(ITCVT1+9 ) + Q(ITCVT2+9 )
               ERR(2,2) = Q(ITCVT1+10) + Q(ITCVT2+10)
               ERR(3,1) = Q(ITCVT1+11) + Q(ITCVT2+11)
               ERR(3,2) = Q(ITCVT1+12) + Q(ITCVT2+12)
               ERR(3,3) = Q(ITCVT1+13) + Q(ITCVT2+13)
               ERR(1,2) = ERR(2,1)
               ERR(1,3) = ERR(3,1)
               ERR(2,3) = ERR(3,2)
               CALL RINV(3,ERR,3,TEMP,IFAIL)
               IF (IFAIL.EQ.0) THEN
                  CHI =
     &                 DX*DX*ERR(1,1)+
     &                 DY*DY*ERR(2,2)+
     &                 DZ*DZ*ERR(3,3)+
     &                 2.0*(
     &                 DX*DY*ERR(1,2)+
     &                 DX*DZ*ERR(1,3)+
     &                 DY*DZ*ERR(2,3))
                  IF (CHI.GT.0.0) THEN
                     P = PROB(CHI,3)
                  ELSE
                     P = 0.0
                  ENDIF

                  IF (P.LT.V3CONF) GOTO 250
               ELSE
                  CALL ERRLOG(IERRL(5),'singular matrix')
                  GOTO 250
               ENDIF
            ENDIF

            NTR = 0
            XGUES(1) = 0.0
            XGUES(2) = 0.0
            XGUES(3) = 0.0

*     add tracks from first vertex
            ITCVT = XXTCVT(IV1)
            DO 920 J=1,TRKSVT(IV1)
               NTR = NTR+1
               TRKS(NTR) = IQ(LQ(ITCVT-1)+1+(J-1)*LENVP)
 920        CONTINUE
            DO 930 K=1,3
               XGUES(K)=XGUES(K)+Q(ITCVT+4+K)/TRKSVT(IV1)
 930        CONTINUE

*     add tracks from second vertex
            ITCVT = XXTCVT(IV2)
            DO 925 J=1,TRKSVT(IV2)
               NTR = NTR+1
               TRKS(NTR) = IQ(LQ(ITCVT-1)+1+(J-1)*LENVP)
 925        CONTINUE
            DO 935 K=1,3
               XGUES(K)=XGUES(K)+Q(ITCVT+4+K)/TRKSVT(IV2)
 935        CONTINUE

*     Create a vertex and momentum bank for this vertex.
            IVRTX = IVRTX + 1
            IF (IVRTX.GT.MAXVRT) THEN
               IVRTX = IVRTX-1
               DO 940 IV=1,IVRTX
                  VON(IV) = 0
 940           CONTINUE
               GOTO 4000
            ENDIF
            CALL MZLIFT(IXTCVX,XXTCVT(IVRTX),XXTCVT(IVRTX),1,MTCVT,0)
*     call desired fitter
            IF(VFITVR.EQ.1) THEN
               CALL TCVRTX(NTR,TRKS,-IVRTX,NDROP,IDROP,
     &              IERR,XGUES)
            ELSEIF (VFITVR.EQ.2) THEN
               CALL TCVRHX(NTR,TRKS,-IVRTX,IERR,XGUES)
            ENDIF

************* CUT3 ************
            IF(IERR.GT.100.OR.NDROP.GT.0) THEN
*     If this was bad, then drop the fit, and forget it ever happened
               CALL MZDROP(IXSTOR,XXTCVT(IVRTX),'L')
               IVRTX=IVRTX-1
               GOTO 250
            ENDIF

*     keep the vertex !!!
*     number of tracksat this vertex
            TRKSVT(IVRTX) = TRKSVT(IV1)+ TRKSVT(IV2)
*     total charge of vertex
            CHARGE(IVRTX) = CHARGE(IV1)+ CHARGE(IV2)
*     total track bit mask of vertex
            MASK(IVRTX)   = MASK(IV1)  + MASK(IV2)
*     turn it on (with flag = 2)
            VON(IVRTX)    = 2
*     number of newly created vertices
            NEWVRT = NEWVRT+1

 250     CONTINUE
 200  CONTINUE

*
* If no new vertices, clean up and exit
*
      IF (NEWVRT.EQ.0) THEN
*     this is the only exit from the routine....
         GOTO 4000
      ENDIF
***************************************************************************
* PART III
***************************************************************************
* Now make valid configuration of vertices
*
* 1) If it can add the vertex to the configuration without overlapping
* with an already present vertex, the current configuration
* is pushed onto the stack.

* 2) If we can push ALL the vertices into a configuration,, the
* configuration is saved.

* 3) Once the vertices run out, the stack is
* popped, and further configurations are tried.

* 4)When the stack is empty, all configurations have been tried
* and the valid ones are given in COLIST()
* this routine should work for any number of vertices or tracks,


* we want the best configuration with lowest vertex multiplicity
* If this doesn't work, we increase MAXNV it until it gets to ntrak
      MAXNV = (OLDVRT+1)/2
+REP,TC_VERT,TCVER3,356.
* start with vertex 1 always
* eventually, we should start from top and work down, it is
* more efficient  FUTURE
      IV = IVRTX
+REP,TC_VERT,TCVER3,362-377.
*     first check if this is a valid vertex
      IF (VON(IV).EQ.0) GOTO 325

*     if the list is full already, then don't both adding anymore
      IF(NVLIST.EQ.MAXNV) GOTO 325

*     also, don't need to go to the very end
      IF(IV+NVLIST.LT.MAXNV) GOTO 325

*     now check to see if we can add it (no overlap)
      IF (IAND(CUMASK,MASK(IV)).EQ.0) THEN
         NVLIST=NVLIST+1
         VLIST(NVLIST) = IV
*     push the state of configuration and last vertex added on the stack
         CALL TCPUSH(CUMASK)
         CALL TCPUSH(IV)
*     this is the current bit mask of tracks used
         CUMASK = IOR(CUMASK,MASK(IV))

*     If the configuration contains all the tracks, the mask will be full
         IF (CUMASK.EQ.FULMSK) THEN
*     jump to end of vertex list
            IV=0
            IF (NCONF.GE.MAXCNF) THEN
               WRITE(CERR,307) NCONF,IVRTX
 307           FORMAT('TCVER3 :Too many vertex configs,Increase MAXCNF',
     &              I5,I5)

               CALL ERRLOG(IERRL(1),CERR)
*     kill vertices, abort, abort, dive, dive, battlestations
               DO 315 IV1=1,IVRTX
                  VON(IV1) = 0
 315           CONTINUE
*     Empty the stack, shut the hatch, raise the mizzenmast.
+REP,TC_VERT,TCVER3,384-400.
*     we don't want it if it has too many particles
               IF(NVLIST.LE.MAXNV) THEN
*     save this configuration
                  NCONF=NCONF+1
                  NVCONF(NCONF) = NVLIST
                  DO 310 IV1=1,NVLIST
                     COLIST(IV1,NCONF) = VLIST(IV1)
 310              CONTINUE
               ENDIF
            ENDIF
         ENDIF
      ENDIF
 325  CONTINUE
*     move down the list of vertices
      IV = IV -1
*     if we are at the bottom, dip into the stack for some more treats
      IF (IV.LE.0) THEN
         CALL TCPOP(IV)
         CALL TCPOP(CUMASK)
*     take that last one off the list
         NVLIST = NVLIST - 1
*     if we are at the bottom of the stack, its time to leave
         IF (CUMASK.EQ.-1) GOTO 350
*     ohoh
         IF (NVLIST.LT.0) THEN
            CALL ERRLOG(IERRL(2),'TCVER3:serious bug NVLIST < 0')
            STOPCB=.TRUE.
            GOTO 4000
         ENDIF
*     go to the next vertex
+REP,TC_VERT,TCVER3,406-665.
* if we got no configuration, lower our sights (raise desired multiplicity)
      IF (NCONF.LE.0) THEN
         IF(MAXNV.LT.(NEWVRT+OLDVRT)) THEN
*     try again, but with higher multiplicity goal
            MAXNV=MAXNV+1
            GOTO 305
         ELSE
            CALL ERRLOG(IERRL(3),
     &        'TCVER3:could not find ANY vertices!,serious bug')
            STOPCB=.TRUE.
            GOTO 4000
         ENDIF
      ENDIF


***************************************************************************
* PART IV
***************************************************************************
*     ok, now we have a list of configurations. We need to test each
*     one and decide which one is the best.

*     find best configuration, based on total chi2
      BESTCH = 1.E6
      BESTCF = 0

*     we now have NV vertices in each configuration
      NV = NVCONF(1)
*     loop over all the configurations
      DO 400 ICONF=1,NCONF
         CHI = 0.0
         NDF = 0

*     first find "primary" vertex (average off all vertices)
         IF (V3BACK.GT.0.0) THEN
            X = 0.0
            Y = 0.0
            Z = 0.0
            DO 430 J=1,NV
               ITCVT1 = XXTCVT(COLIST(J,ICONF))
               X = X+Q(ITCVT1+5)
               Y = Y+Q(ITCVT1+6)
               Z = Z+Q(ITCVT1+7)
 430        CONTINUE
            X = X/NV
            Y = Y/NV
            Z = Z/NV
         ENDIF

*     loop over all the vertices in this configuration
         DO 410 J=1,NV
            ITCVT1 = XXTCVT(COLIST(J,ICONF))
            CHI = CHI +  Q(ITCVT1+14)
            NDF = NDF + IQ(ITCVT1+4)

*     optional step: if momentum of vertex is pointing towards
*     primary vertex, then add a penalty of V3BACK to chi
*     dx,dy,dz is displacement from primary vertex
*     px,py,pz is total momentum at the vertex
            IF (V3BACK.GT.0.0) THEN
               DX = Q(ITCVT1+5)-X
               DY = Q(ITCVT1+6)-Y
               DZ = Q(ITCVT1+7)-Z
*     calculate total momentum at vertex
               PX = 0.0
               PY = 0.0
               PZ = 0.0
               ITCVP = LQ(ITCVT1-1)
               DO 420 K=1,IQ(ITCVT1 + 1 )
                  PX =  PX + Q(ITCVP+4)
                  PY =  PY + Q(ITCVP+5)
                  PZ =  PZ + Q(ITCVP+6)
                  ITCVP = ITCVP + LENVP
 420           CONTINUE
               IF (DX*PX+DY*PY+DZ*PZ.LT.0.0) CHI=CHI+V3BACK
            ENDIF

 410     CONTINUE

*     protect from divide by zeros
         IF (NDF.LE.0) NDF = 1

*     calculate reduced chi^2
         CHI = CHI/NDF

*     save this configuration if it is the best
         IF (CHI.LT.BESTCH) THEN
            BESTCF = ICONF
            BESTCH = CHI
         ENDIF
 400  CONTINUE

*     turn off all vertices, for now, but we turn them on again soon
      DO 450 IV=1,IVRTX
         VON(IV) = 0
 450  CONTINUE

*     check to see if we have a configuration
      IF (BESTCF.EQ.0) THEN
         CALL ERRLOG(IERRL(4),'TCVER3:No valid configs==0.0')
*     forget this event, clean up...
         STOPCB=.TRUE.
         GOTO 4000
      ENDIF

*     This section can be used for debugging.
*     Uncomment the desired BESTCF = ... line
C*     force the first one
C      IF(NCONF.EQ.2) THEN
C*     worst
C         BESTCF = 3 - BESTCF
C*     first
C         BESTCF = 1
C*     second
C         BESTCF = 2
C*     best (do nothing)
C         BESTCF = BESTCF
C      ELSE
C* drop event
C         GOTO 4000
C      ENDIF

*     turn on selected vertices
      DO 460 IV=1,NV
         VON(COLIST(IV,BESTCF)) = 1
 460  CONTINUE
      OLDVRT = NV

*     try again for smaller multiplicity
*     For debuggin, comment this line out, to force only 2-preong vertices
      GOTO 1000

***************************************************************************
* PART V
***************************************************************************

*     clean up, finish everything
+REP,TC_VERT,TCVER3,671-679.
      DO 470 IV=1,IVRTX
         IF (VON(IV).EQ.0) THEN
            CALL MZDROP(IXSTOR,XXTCVT(IV),'L')
         ELSE
*     good vertex, move the bank to the correct spot
            IGOOD=IGOOD+1
            CALL ZSHUNT(IXSTOR,XXTCVT(IV),LTCVX,-IGOOD,0)
            ITCVT = XXTCVT(IV)

*     assign an OFFICIAL vertex # (IGOOD) to all the tracks

+REP,TC_VERT,TCVER3,702-725.
*     that's it.


***************************************************************************
* PART VI
***************************************************************************
*     now we just make sure that the bank is totally corrupt due to my
*     stupidity
9999  CONTINUE

      NV=IQ(LTCVX+1)

      IF(NV.NE.IVRTX) THEN
         CALL ERRLOG(IERRL(6),'TCVER3:vertex # check')
         STOPCB=.TRUE.
      ENDIF

      IF(NV.GT.0) THEN
         DO IV=1,NV
            IF(LQ(LTCVX-IV).LE.0) THEN
               CALL ERRLOG(IERRL(7),'empty link to TCVT')
               STOPCB=.TRUE.
            ELSE
               ITCVT1=LQ(LTCVX-IV)
               ITCVP=LQ(ITCVT1-1)
               IF(ITCVP.LE.0) THEN
                  CALL ERRLOG(IERRL(8),'empty link to TCVP')
                  STOPCB=.TRUE.
               ELSE
                  ITCVH=LQ(ITCVP)
                  IF(ITCVH.LE.0) THEN
                     CALL ERRLOG(IERRL(9),'empty link to TCVH')
                     STOPCB=.TRUE.
                  ENDIF
               ENDIF
            ENDIF
         ENDDO
+ADD,TC_VERT,TCVRHX,0.
*CMZU:  2.01/13 14/03/97  16.19.51  by  M.T. Lakata
*CMZU:  2.01/11 29/01/97  17.41.31  by  M.T. Lakata
*CMZU:  2.01/10 22/01/97  15.22.16  by  M.T. Lakata
*CMZU:  2.01/09 26/11/96  19.49.18  by  Rafik Ouared
*CMZ :          02/10/96  17.27.42  by  Rafik Ouared
*               only acceptable momenta kept.
*CMZU:  2.01/07 17/06/96  17.28.46  by  F.-H.Heinsius
*CMZU:  2.01/05 02/04/96  15.41.22  by  F.-H.Heinsius
*CMZ :  2.01/01 13/03/96  15.30.55  by  F.-H.Heinsius
+REP,TC_VERT,TCVRHX,14-26.
*     Modified:      27 March 1996  MTL
*                    Added multiple scattering, fixed many bugs,
*                    added outlier rejection, fixed momentum errors,
*                    etc,etc,etc
*     Modified:      02 October 1996
*                    Job was crashing because very large momenta in SPP matrix.
*                    max value for SPP(i,j) set to 1.D+10.
*
*                    14 March 1997 MTL
*                    (1) zebra banking lift was sloppy, now fixed to allow for
*                    garbage collection.
*                    (2) sigma(r) removed as input to the fit because there
*                    is essentially no constraint on this value. Tests show
*                    that there is no significant change in the efficiency,
*                    while the chi^2 makes more sense now, since the original
*                    position measurement is based on only on one variable
*                    (time).
*
*
*     References:    Siegmund Brandt, Data Analysis, Chapter 9
*     editions: 1970, 1976 or 1992
*
*     Description: This routine will attemps to refit the tracks
*     listed in NTRAK, with the extra constraint that they
*     all pass through a common vertex.  It does a similar
*     job as TCVRTX, returning the same zebra banks and
*     error codes(and 900).  However, instead of simply fitting
*     a vertex to the helices, it refits the helices
*     altogether, so the routine is more similar to TCHELX,
*     except instead of 5*NTRKS free parameters for the
*     helices, it uses only 3+3*NTRKS.
+REP,TC_VERT,TCVRHX,33-41.
*     and:    _Datenanalyse or Data Analysis_, Siegmund Brandt
*             p. 183  (1st ed.)
*             p. 271? (2nd ed.)
*             p. 318  (3rd ed.)
*
*     Called by:     TCVER3
*
*     Subroutines Referenced:
*     DINV (cernlib F010).
*     TXTGET (locater) gets hits
*     MZLIFT (cernlib, zebra)
*     ERRLOG (cboff) error message logging
+REP,TC_VERT,TCVRHX,66-116.
*     Word    Type      Description
*     ================================================
*     LQ(-1)   I*4       Pointer to the TCVP data bank
*     ------------------------------------------------
*     +1      I*4       Number of tracks to this vertex.
*     +2      I*4       Number of neutrals.
*     +3      I*4       Error code of fit.
*     +4      I*4       Number of degrees of freedom.
*     +5      R*4       x-vertex  (cm).
*     +6      R*4       y-vertex  (cm).
*     +7      R*4       z-vertex  (cm).
*     +8      R*4       covarianve matrix for vertex.
*
*     +14     R*4       Chisquare of the fit.
*
*     The format of the TCVP bank is:
*
*     Word    Type      Description
*     -------------------------------------------------
*     +1      I*4       Track number in TCTR
*     +2      I*4       Quality word.
*     +3      R*4       Charge of track.
*     +4      R*4       x-momentum (MeV/c)
*     +5      R*4       y-momentum (MeV/c)
*     +6      R*4       z-momentum (MeV/c)
*     +7      R*4       Energy (assuming pion mass) (Mev).
*     +8      R*4       Covariance matrix for the momentum
*     +9      R*4          "
*     +10     R*4          "   |  8       |   | s_xx           |
*     +11     R*4          "   |  9 10    | = | s_xy s_yy      |
*     +12     R*4          "   | 11 12 13 |   | s_xz s_yz s_zz |
*     +13     R*4          "
*     === Repeat the last 13 for every track at this vertex ===
*
*
*     Bank Name          : TCVH    "T"ra"C"king,"V"ertex fitted "H"elix
*     PREVIOUS bank link : ITCVP    i.e. ITCVH = LQ(ITCVP)
*
*     TCVH bank format:
*
*     Word  Type  Description    equivalent to TCTX word
*     -------------------------------------------
*     +1    R*4  r0                +11  --
*     +2    R*4  z0                +12   |
*     +3    R*4  alpha             +13   |
*     +4    R*4  tan lambda        +14   |- first track at vertex
*     +5    R*4  psi0              +15   |
*     +6    R*4  s                 +16   |
*     +7    R*4  Chi**2            +17  --
*     The structure then repeats every 7=LENVH words for the rest of the
*     tracks at that vertex. i.e. Offset +8 is the r0 of track 2.
+REP,TC_VERT,TCVRHX,131-140.
&SEQ,TCV3CM.

C+SEQ,TCSCAT. (for GASXTC)

*
*     CONSTANTS: --------------------------------------------------------
*
*     MAXTRK  maximum number of tracks
*     MAXPT   maximum number of hits on a track. Should be same as TXTGET
*     MAX2PT  MAXPT*2
*     MAX3PT  MAXPT*3
*
      INTEGER MAXTRK,MAXPT,MAX2PT,MAX3PT,MAXPAR
      PARAMETER (MAXTRK=10,MAXPT=50,MAX2PT=MAXPT*2,MAX3PT=MAXPT*3)
+REP,TC_VERT,TCVRHX,158.
*     Multiple scattering:
*     SIGRF Is an assumed measurement error of 130 microns squared times 720.
*     SIGZ  Is an assumed measurement error of 7mm squared times 12.
*     GASXT is inverse scattering length in JDC gas
*        Note: GASXT is recalculated by Mark Lakata (it does not include wires)
*              In principle GASXTC (Common TCSCAT)should be identical
*
      REAL   SIGRF,SIGZ,GASXT
      PARAMETER (SIGRF=720.*0.0130**2,SIGZ=12.0*0.7**2,GASXT=5.5E-5)

*     for "1-prong vertices", chose some reasonable errors
*     ERLONG = longitudinal error (also z error?), 1 cm
*     ERLAT  = transverse error, 1 mm
      DOUBLE PRECISION ERLONG, ERTRAN
      PARAMETER (ERLONG = 1.D0, ERTRAN = 1.D-1)
+REP,TC_VERT,TCVRHX,170-174.
*     IERR1  returned error code.
*     XGUES is a guess to the vertex psoition.
*
      INTEGER  NTRKS,NTRAK(*),ITCVT1,IERR1
+REP,TC_VERT,TCVRHX,213.
      INTEGER     I,J,K,J1,J2,KX,KY,KZ,IX,IY,IZ,II,JJ,KK
+REP,TC_VERT,TCVRHX,256-262.
      DOUBLE PRECISION X(MAXPAR),XSI(MAXPAR),X0(MAXPAR)
      DOUBLE PRECISION ATGBA(MAXPAR,MAXPAR),ATGBC(MAXPAR)
      DOUBLE PRECISION Y(MAX3PT,MAXTRK),Y0(MAX3PT,MAXTRK)
      DOUBLE PRECISION GYINV(MAX3PT,MAXTRK),C(MAX2PT,MAXTRK)
      DOUBLE PRECISION GB(3,MAXPT,MAXTRK),ATGB1,ATGB2
      DOUBLE PRECISION AJ1(6),AJ2(6)
      DOUBLE PRECISION A(MAX2PT,6,MAXTRK),B(MAX2PT,3,MAXTRK)
      DOUBLE PRECISION DC1,DC2,DY1,DY2,DY3
+REP,TC_VERT,TCVRHX,283-286.
*     Error message numbers needed for cboff routine ERRLOG
      INTEGER IERRL(10)
      CHARACTER*160 CERR

*     RR        = radial position of hit, from center of detector
*     P         = phi position of hit, from center of detector
*     COSP,SINP = cos,sin of Phi
*     DYP       = update to P position
*     NPTS0     = points from txtget, .LE. NPTS(IT)

      DOUBLE PRECISION RR,COSP,SINP,DYP
      INTEGER NPTS0

*     CH1,CH2,CHI12= chi^2 matrix of hits for two constraints
*     BADPNT    = If True, this point should be rejected
*     BADCOUNT  = number of badpoints on this track
      DOUBLE PRECISION CH2,CH1,CH12
      LOGICAL BADPNT(MAXPT,MAXTRK)
      INTEGER BADCOUNT(MAXTRK)

*     This is kept in case one wants to do checks on the chi^2 of the
*     individual hits
*     The first parts can be put in a USINIT file, to generate
*     a histogram
c$$$#ifdef MTLCHI
c$$$      character*8 addch(15)
c$$$      data addch/ 'rr','p','it','chp','chpr','chr','ru','ev','ipt',
c$$$     &     'dyp','dyr','sp','sr','dyz','sz'/
c$$$
c$$$      CALL HLIMIT(-pawsize)
c$$$      CALL HROPEN(1,'VERTEX',filename,'NCXP',1024,ISTAT)
c$$$      call hbookn(4,'addchi',15,'//VERTEX',8192,addch)
c$$$      write(*,*) 'hbookn 4 done'
c$$$#endif

*     this part and the part later in the file stay here.
c$$$#ifdef MTLCHI
c$$$*     P   = phi of hit
c$$$      DOUBLE PRECISION P
c$$$* only for testing purposes...
c$$$      real nt(100)
C$$$      integer nmtlchi /0/
C$$$      save nmtlchi
c$$$#endif
*        LYR0  First layer of the track.
*        LYRN  Last  layer of the track.
*
      INTEGER  LYR0,LYRN,NPTSX
*        The following 4 parameters come from Gluckstern.
*
      REAL    CN,DN,EN,FN
*
      PARAMETER (CN = 1.43)
      PARAMETER (DN = 0.21)
      PARAMETER (EN = 0.23)
      PARAMETER (FN = 0.20)
*
*        XLEN   Distance in [cm] that the particle travels in the JDC gas.
*        ALM    is (14.1/p*beta)**2 (beta = v/c).
*
      REAL    XLEN,ALM
*
*        SECL  is the sec(lambda)
*        SECL2 is sec**2(lambda)
*        PINV  is 1/p(total)**2
*        PTOT  is p(total)**2
*        ETOT  is E**2   (assuming a pion).
*
      REAL    SECL,SECL2,PINV,PTOT,ETOT
*
*        The following S values are the corrections to the covariance matrix.
*        These are equivalenced to a 15 element array, and added to the stored
*        matrix in a do loop.
*
      REAL    SAA,S11,S22,S33,S44,S35,S55

      REAL R0X(6),R2X(6),SCTRAD(0:6),SCTLEN(6),ALMGAS
      INTEGER IRHO,IRHO2
      SAVE R0X,R2X,SCTRAD
+REP,TC_VERT,TCVRHX,296.
*     this may be uncommented later, if the updated hits are desired
+ADD,TC_VERT,TCVRHX,308.
*     These are the radiation lengths and radii of the various
*     scattering materials between r=0 and the first layer of the JDC
*
*                       H2     Mylar   Al      PWC1    PWC2    JDC
      DATA SCTLEN /    1.0E-3, 0.3E-3, 5.6E-3, 1.1E-3, 1.1E-3, 2.5E-3/

*     Radii in [cm]
*                       H2     Mylar   Al      PWC1    PWC2    JDC
      DATA SCTRAD /0.0, 0.0,   0.85,   1.7,    2.6,    4.3,    5.9/
+REP,TC_VERT,TCVRHX,315-354.
*     If any changes to defaults, print out warning, otherwise
*     WELCOME to TCVRHX !

         WRITE(LLOG,294)
 294     FORMAT(' *** TCVRHX Initialized, Version 1.1')

*     check the other parameters
         CALL TCVPCR('VHMCH',VHMCH,VHMCHD)
         CALL TCVPCR('VHGCH',VHGCH,VHGCHD)
         CALL TCVPCR('VHCFL',VHCFL,VHCFLD)
         CALL TCVPCR('VHERC',VHERC,VHERCD)
         CALL TCVPCI('VHBDT',VHBDT,VHBDD)
         CALL TCVPCR('VHORT',VHORT,VHORD)
         CALL TCVPCR('VHSPT',VHSPT,VHSPD)
         CALL TCVPCR('VHSZT',VHSZT,VHSZD)

*     Setup for multiple scattering
*     1)correct the scattering length with the log correction
*     2)calculate the sum X_i        for each scatterer j (i>=j) (R0X)
*     3)calculate the sum X_i D_ij^2 for each scatterer j (i>=j) (R2X)
*
*     X_i is the radiation length of the scatter
*     D_ij is the distance between scatterers i and j

         DO 50 IRHO=6,1,-1
            SCTLEN(IRHO) = SCTLEN(IRHO)*(1.0+0.038*LOG(SCTLEN(IRHO)))
            IF (IRHO.EQ.6) THEN
               R0X(IRHO) = SCTLEN(IRHO)
            ELSE
               R0X(IRHO) = R0X(IRHO+1)+SCTLEN(IRHO)
            ENDIF
 50      CONTINUE

         DO 51 IRHO=1,6
            R2X(IRHO) = 0.0
            DO 52 IRHO2=IRHO,6
               R2X(IRHO) = R2X(IRHO)+
     &              SCTLEN(IRHO2)*(SCTRAD(IRHO2)-SCTRAD(IRHO-1))**2
 52         CONTINUE
 51      CONTINUE

      ENDIF
*     set status bits in TCVX header word
*     bits 1-8  tcver[123]   :
*     bits 9-16 tcvr[th]x    : 0 = tcvrtx
*                              1 = tcvrhx
      CALL SBYT(1,IQ(LTCVX+0),9,8)

*     If ITCVT1 is negative, then it is the index of XXTCVT common
*     if it is positive, then it is an existing vertex number

      IF (ITCVT1.LT.0) THEN
         ITCVT = XXTCVT(-ITCVT1)
      ELSE
         ITCVT = LQ(LTCVX-ITCVT1)
+REP,TC_VERT,TCVRHX,364.
 296     FORMAT('<TCVRHX> called with ',I3,' tracks, out of range ',
+ADD,TC_VERT,TCVRHX,395.
         RCHISQ = 0.0D0
+DEL,TC_VERT,TCVRHX,437-451.
+DEL,TC_VERT,TCVRHX,468-475.
+REP,TC_VERT,TCVRHX,488.

*     Now for the error/matrix.  We know that the "vertex" is on the
*     track, but we don't know where. Thus we put a long skinny error
*     ellipse on the best guessed point.  If the track ends in the jdc,
*     it probably began between layers, so let's assign about 1 cm to
*     the position. If it ends inside the pwc, then it probably was
*     inside the target too, so let's also give it 1 cm errors.  These
*     errors are only to keep things from being sane, and are not meant
*     to be accurate, since there is absolutely no way of knowing where
*     the vertex is without another track's information.
*
*     ERTRAN = transverse error (typically 1 mm)
*     ERLONG = longitudinal error (typically 1 cm)

*     We will use sensible errors, and one that is stretched along the
*     track.  First we make an ellipse with major axis on x-axis, and
*     rotate the axes so that the major axis lies along the track
*     direction (angle P).  This requires transforming the error matrix
*     by the usually tensor transformation.
*
*     axes rotation:
*     x' = cosp x  -sinp y
*     y' = sinp x  +cosp y
*
*     tensor transformation: (note 12 term is zero,so it is not included)
*     11  = cos cos 11 + sin sin 22
*     12  = cos sin 11 +-sin cos 22
*     22  = sin sin 11 + cos cos 22
*
*     we finally use the fact that p = beta +- 90s, but p = p +- 180,
*     so we can use the following relations:
*
         COSP = SINB
         SINP = -COSB
         ATGBA(1,1) = (COSP*ERLONG)**2 + (SINP*ERTRAN)**2
         ATGBA(1,2) = COSP*SINP*(ERLONG**2-ERTRAN**2)
         ATGBA(2,1) = ATGBA(1,2)
         ATGBA(2,2) = (SINP*ERLONG)**2 + (COSP*ERTRAN)**2
         ATGBA(3,3) = ERLONG**2

*     Fill vertex bank with result
+REP,TC_VERT,TCVRHX,513-555.
         CALL TXTGET(12,NTRAK(IT),NPTS0,R1,X1,Y1,Z1,D1,D2,D3,
     &        ILAY,ISC,IRS,TIM,AMPL,AMPR,SR,SP)

         NPTS(IT) = 0
         BADCOUNT(IT)= 0
         IX = 1
         DO 505 IPT=1,NPTS0
            IY = IX +1
            IZ = IX +2

*     Reject the outliers COMPLETELY, don't just assign a big error to
*     them and hope they go away, since this messes up the confidence
*     level at the end.  Except for PWC z-errors, since this are
*     unknown....
*     Outlier z errors are pegged at 86 cm. remove these points.
*     pwc errors are 21.5 cm.. be careful with these.
*
            IF(SP(IPT).LT.0.6.AND.(IPT.EQ.1.OR.D3(IPT).LT.20.0)
     &           .AND.D3(IPT).LT.80.0) THEN
*     MTL: scale errors by some amount
               SP(IPT) = SP(IPT)*VHSPT
               D3(IPT) = D3(IPT)*VHSZT
*     Put lower bound on radius error of 1 micron (approx width of sense
*     wire)

*     phi
               Y0   (IX,IT) = DBLE(X1(IPT))
               Y    (IX,IT) = Y0(IX,IT)
               GYINV(IX,IT) = DBLE(VHERC*SP(IPT)**2)
*     radius
               Y0   (IY,IT) = DBLE(Y1(IPT))
               Y    (IY,IT) = Y0(IY,IT)
* error on radius is unused
               GYINV(IY,IT) = 0.0
*     zed
               Y0   (IZ,IT) = DBLE(Z1(IPT))
               Y    (IZ,IT) = Y0(IZ,IT)
               GYINV(IZ,IT) = DBLE(VHERC*D3(IPT)**2)
               NPTS(IT) = NPTS(IT) +  1
               IX = IX + 3
               BADPNT(NPTS(IT),IT) = .FALSE.
            ENDIF
 505     CONTINUE

         IF (NPTS(IT).LT.V3MNHT.OR.NPTS(IT).EQ.0) THEN
            IERR1  = 900
C     WRITE(*,*) 'TCVRHX: No good points on this track ',
C     &              NPTS0,NPTS(IT),IEHDCB(4),IEHDCB(5)
            GOTO 5000
         ENDIF
         TOTPT = TOTPT + NPTS(IT)
 200  CONTINUE

*     The Number of degrees of freedom in this fit is
*     N_df = N_c * N_pt - N_par
*     N_c  = number of constraint equations = 2
*     N_pt = number of hits
*     N_par= number of fit parameters
*
         NPAR  = 3 + 3*NTRKS
         NDF   = 2*TOTPT - NPAR

*     Now we need to have an initial guess for the common vertex.
*     Since TCVER3 will be handing us 2-prong vertices, this will be
*     a good start...

         IF (NTRKS.EQ.2) THEN
+REP,TC_VERT,TCVRHX,563-571.
*     The equation of a line that passes through (x,y)a and (x,y)b
*     is
*     (yb-ya) x  + (xa-xb) y = (xa yb - ya xb)
*
*     Thus we need to solve the system, which is the
*     intersection of two tracks
*
*     (yb1-ya1) x  + (xa1-xb1) y = (xa1 yb1 - ya1 xb1)
*     (yb2-ya2) x  + (xa2-xb2) y = (xa2 yb2 - ya2 xb2)
+REP,TC_VERT,TCVRHX,577-672.
            DO 4564 IT =1,2
               AA(IT,1) = Y(YB,IT) - Y(YA,IT)
               AA(IT,2) = Y(XA,IT) - Y(XB,IT)
               BB(IT)   = Y(XA,IT)*Y(YB,IT)-Y(YA,IT)*Y(XB,IT)
 4564       CONTINUE

*     invert AA

            DET = AA(1,1)*AA(2,2)-AA(1,2)*AA(2,1)
            IF(DET .EQ. 0.D0) THEN

*     An attempt was made to invert a singular matrix.

               CALL ERRLOG(IERRL(4),'<TCVRHX> AA is singular')
               IERR1 = 700
               GOTO 5000
            ENDIF

            AAINV(1,1) =  AA(2,2)/DET
            AAINV(2,2) =  AA(1,1)/DET
            AAINV(1,2) = -AA(1,2)/DET
            AAINV(2,1) = -AA(2,1)/DET
*
            X(VX) = DBLE(AAINV(1,1)*BB(1)+AAINV(1,2)*BB(2))
            X(VY) = DBLE(AAINV(2,1)*BB(1)+AAINV(2,2)*BB(2))
*     now that we know vx, we can extrapolate the Z position at vx
*     from each track, and then just average the two together. The z position
*     is not as important as the (x,y) position.

            IF (AA(1,2).NE.0D0.AND.AA(2,2).NE.0D0) THEN
               X(VZ) = (
     &              Y(ZA,1)-(Y(ZB,1)-Y(ZA,1))*(X(VX)-Y(XA,1))/AA(1,2)+
     &              Y(ZA,2)-(Y(ZB,2)-Y(ZA,2))*(X(VX)-Y(XA,2))/AA(2,2)
     &              )*0.5D0
            ELSEIF (AA(1,1).NE.0D0.AND.AA(2,1).NE.0D0) THEN
               X(VZ) = (
     &              Y(ZA,1)-(Y(ZB,1)-Y(ZA,1))*(X(VY)-Y(YA,1))/AA(1,1)+
     &              Y(ZA,2)-(Y(ZB,2)-Y(ZA,2))*(X(VY)-Y(YA,2))/AA(2,1)
     &              )*0.5D0
            ELSE
               X(VZ) = 0.D0
            ENDIF
            IF (X(1)**2+X(2)**2.GT.2.D0.AND.DABS(X(3)).GT.2.D0) THEN
               DO 111 I=1,3
                  IF (DABS(X(I)).GT.30.D0.OR.
     &                 DABS(X(I)).GT.DABS(Y(I,1)).OR.
     &                 DABS(X(I)).GT.DABS(Y(I,2))) THEN
*
*     The guess was bad, but don't give up... try 0,0,0 !
*
                     X(1) = 0.D0
                     X(2) = 0.D0
                     X(3) = 0.D0
                  ENDIF
 111           CONTINUE
            ENDIF
         ELSE
            X(VX) = DBLE(XGUES(1))
            X(VY) = DBLE(XGUES(2))
            X(VZ) = DBLE(XGUES(3))
         ENDIF

*     Initialization.
*     set Chi-square to a very large value.
*
         OLDCH = 1.0D10
         DELCH = 1.0D10
         ITER  = 1

*     Load the helix fit information from the TCTR bank for each
*     track.
*
         IOFF = 0
         DO 800 IT = 1,NTRKS
            ITCTR = LQ(LTCTR-NTRAK(IT))
            IOFF  = IOFF  +  3
*
*     initial guesses
            SGN(IT)        = DBLE(Q(ITCTR+16))

*     alpha, parameter1
            X(ALPHAX+IOFF) = DBLE(Q(ITCTR+13))

*     tan (lambda) parameter2
            X(TANLAX+IOFF) = DBLE(Q(ITCTR+14))

*     phi0, parameter3
            RP = (SGN(IT)*DBLE(Q(ITCTR+11))+1/X(ALPHAX+IOFF))
            PHI0 = DATAN2(
     &           X(2)-RP*DSIN(DBLE(Q(ITCTR+15))-SGN(IT)*DPIHLF),
     &           X(1)-RP*DCOS(DBLE(Q(ITCTR+15))-SGN(IT)*DPIHLF))
            IF (PHI0.LT.0.0) PHI0 = PHI0 + DPI2
            X(PHI0XX+IOFF) = PHI0
*
 800     CONTINUE



*     Save the initial values (for calculation of chi^2)

         DO 905 I=1,NPAR
            X0(I) = X(I)
 905     CONTINUE
+REP,TC_VERT,TCVRHX,679-715.
*     beginning of main iteration loop.
*     control comes back here if the first pass is not good enough

 1000    CONTINUE

*     zero some things

         CHISQ  = 0.0D0
         DO 1100 J = 1,NPAR
            XSI(J)   = 0.D0
            ATGBC(J) = 0.D0
            DO 1120 K=1,NPAR
               ATGBA(J,K) = 0.D0
 1120       CONTINUE
 1100    CONTINUE
         IOFF = 0

         DO 2800 IT = 1,NTRKS

*     Compute repeatedly used quantities in the loop:

            IOFF    = IOFF + 3

            ALPHA   = X(ALPHAX+IOFF)
            TANLA   = X(TANLAX+IOFF)
            PHI0    = X(PHI0XX+IOFF)
            SGN1    = SGN(IT)
            IF(ALPHA.NE.0D0) THEN
               INVALP  = 1.D0/ALPHA
            ELSE
               INVALP  = 1.D20
            ENDIF
            COSPHI  = DCOS(PHI0)
            SINPHI  = DSIN(PHI0)
+ADD,TC_VERT,TCVRHX,722.
            IF(BADPNT(IPT,IT)) GOTO 2000
+REP,TC_VERT,TCVRHX,756-772.
            IF(BETA.LT.0.0D0.AND.ITER.GT.1) THEN
C               IF (IPT.GT.1) THEN
C                  IF (BETA-OLDBET.LT.0.D0) THEN
               WRITE(CERR,309) SNGL(BETA),SNGL(PHI0),
     &              NINT(SNGL(SGN1)),
     &              IT,IPT,SNGL(DATAN2(SINB,COSB))
 309           FORMAT('<TCVRHX> Beta <0',
     &              'Beta=',F10.3,'phi',F10.3,'sign',I2,
     &              'track',I2,'hit',I2,'atan',F10.3)
               CALL ERRLOG(IERRL(3),CERR)

               IERR1 = 700
               GOTO 5000

            ENDIF
C            ENDIF
+REP,TC_VERT,TCVRHX,785-798.
            RR   = SQRT(XX*XX+YY*YY)
            IF(RR.NE.0D0) THEN
               COSP = XX/RR
               SINP = YY/RR
            ELSE
               COSP = 1.0D0
               SINP = 0.0D0
            ENDIF
*     dF1/dp =
*     dF1/dr =
*     dF1/dz = 0
            B(J1,1,IT) = 2.D0 * RR*(-COSB*SINP+SINB*COSP)
C            B(J1,2,IT) = 2.D0 *    ( COSB*COSP+SINB*SINP)
            B(J1,3,IT) = 0.D0

*     dF2/dp =
*     dF2/dr =
*     dF2/dz = -1.0
            B(J2,1,IT) = -SGN1*TANLA*RR*(SINB*SINP+COSB*COSP)/C2S2B
C            B(J2,2,IT) =  SGN1*TANLA   *(SINB*COSP-COSB*SINP)/C2S2B
+REP,TC_VERT,TCVRHX,807-812.
C            GB11 = GYINV(KX,IT)*B(J1,1,IT)*B(J1,1,IT)
C     &           + GYINV(KY,IT)*B(J1,2,IT)*B(J1,2,IT)
C            GB12 = GYINV(KX,IT)*B(J1,1,IT)*B(J2,1,IT)
C     &           + GYINV(KY,IT)*B(J1,2,IT)*B(J2,2,IT)
C            GB22 = GYINV(KX,IT)*B(J2,1,IT)*B(J2,1,IT)
C     &           + GYINV(KY,IT)*B(J2,2,IT)*B(J2,2,IT)
C     &           + GYINV(KZ,IT)
            GB11 = GYINV(KX,IT)*B(J1,1,IT)*B(J1,1,IT)
            GB12 = GYINV(KX,IT)*B(J1,1,IT)*B(J2,1,IT)
            GB22 = GYINV(KX,IT)*B(J2,1,IT)*B(J2,1,IT)
+REP,TC_VERT,TCVRHX,850-851.
            A(J1,VX,      IT) = -2.D0 * COSB
            A(J1,VY,      IT) = -2.D0 * SINB
+REP,TC_VERT,TCVRHX,859-860.
            A(J2,VX,      IT) = -SGN1 * TANLA * SINB/C2S2B
            A(J2,VY,      IT) =  SGN1 * TANLA * COSB/C2S2B
+REP,TC_VERT,TCVRHX,868.
CMARK vrhx3.4

            DO 1860 K=1,6
               AJ1(K) = A(J1,K,IT)
               AJ2(K) = A(J2,K,IT)
 1860       CONTINUE
+REP,TC_VERT,TCVRHX,878-879.
               ATGB1 = AJ1(K)*GB11 + AJ2(K)*GB12
               ATGB2 = AJ1(K)*GB12 + AJ2(K)*GB22
+REP,TC_VERT,TCVRHX,890.
     &                 + ATGB1*AJ1(J) + ATGB2*AJ2(J)
+REP,TC_VERT,TCVRHX,896-897.
     &              + ATGB1 * C(J1,IT) + ATGB2*C(J2,IT)
 1800       CONTINUE
CMARK vrhx3.44
+REP,TC_VERT,TCVRHX,937-938.
            I = ALPHAX+IOFF
            IF (DSQRT(ATGBA(I,I)).GT.DABS(X(I))) THEN
+ADD,TC_VERT,TCVRHX,950.
            IF(BADPNT(IPT,IT)) GOTO 2700
+REP,TC_VERT,TCVRHX,975-978.
            DYP  = GYINV(KX,IT) *(B(J1,1,IT) *FACT0 + B(J2,1,IT)*FACT1)
C            DYR  = GYINV(KY,IT) *(B(J1,2,IT) *FACT0 + B(J2,2,IT)*FACT1)
            XX   = Y(KX,IT)
            YY   = Y(KY,IT)
            RR   = SQRT(XX*XX+YY*YY)
            IF (RR.NE.0.D0) THEN
               COSP = XX/RR
               SINP = YY/RR
            ELSE
               COSP = 1.0D0
               SINP = 0.0D0
            ENDIF
C            DY1  = -RR*SINP*DYP + COSP*DYR
C            DY2  =  RR*COSP*DYP + SINP*DYR
            DY1  = -RR*SINP*DYP
            DY2  =  RR*COSP*DYP
            DY3  = GYINV(KZ,IT) *(0.D0              - 1.0D0     *FACT1)
+REP,TC_VERT,TCVRHX,996-1005.
*     inefficient, but necessary for the tests. perhaps this can be
*     fixed.

*     change in phi, since start of fit
            DYP = DATAN2(Y (KY,IT),Y (KX,IT)) -
     &            DATAN2(Y0(KY,IT),Y0(KX,IT))
*     change in radius, since start of fit
C            DYR = SQRT(Y (KX,IT)**2+Y (KY,IT)**2)-
C     &            SQRT(Y0(KX,IT)**2+Y0(KY,IT)**2)
C            FACT0 = B(J1,1,IT)*DYP + B(J1,2,IT)*DYR
C            FACT1 = B(J2,1,IT)*DYP + B(J2,2,IT)*DYR - DY3
            FACT0 = B(J1,1,IT)*DYP
            FACT1 = B(J2,1,IT)*DYP - DY3
*
*
            CH1  = GB(1,IPT,IT)*FACT0*FACT0
            CH12 = GB(2,IPT,IT)*FACT0*FACT1
            CH2  = GB(3,IPT,IT)*FACT1*FACT1

            IF(ITER.EQ.2.AND.VHBDT.GT.0) THEN
*     LET's try throwing out outliers, based on ch1, only on second
*     pass. The first point can not be thrown out.
               IF (CH1.GT.VHORT.AND.IPT.GT.1) THEN

                  BADPNT(IPT,IT) = .TRUE.
*     count the number of bad hits on this track
                  BADCOUNT(IT) = BADCOUNT(IT)+1
*     reduce the Ndf (2 = 2 contraint eqns)
                  NDF = NDF - 2
                  CH1  = 0.0D0
                  CH12 = 0.0D0
                  CH2  = 0.0D0
                  IF (NPTS(IT)-BADCOUNT(IT).LT.V3MNHT.OR.
     &                 NDF.LT.1) THEN
                     IERR1 = 400
                     GOTO 5000
                  ENDIF
               ENDIF
            ENDIF
c$$$#ifdef MTLCHI
c$$$* only histogram on the second iteration
c$$$            if(iter.eq.3.and.nmtlchi.lt.100000) then
c$$$               nmtlchi = nmtlchi + 1
c$$$               rr = SQRT(Y0(KX,IT)**2+Y0(KY,IT)**2)
c$$$               p  = DATAN2(Y0(KY,IT),Y0(KX,IT))
c$$$
c$$$               nt(1) = sngl(rr)
c$$$               nt(2) = sngl(p)
c$$$               nt(3) = float(it)
c$$$               nt(4) = sngl(ch1)
c$$$               nt(5) = sngl(ch12)
c$$$               nt(6) = sngl(ch2)
c$$$               nt(7) = float(iehdcb(4))
c$$$               nt(8) = float(iehdcb(5))
c$$$               nt(9) = float(ipt)
c$$$               nt(10)= sngl(dyp)
c$$$C               nt(11)= sngl(dyr)
c$$$               nt(12)= sngl(GYINV(KX,IT))
c$$$C               nt(13)= sngl(GYINV(KY,IT))
c$$$
c$$$* z pulls
c$$$               nt(14)= sngl(dy3)
c$$$               nt(15)= sngl(GYINV(KZ,IT))
c$$$               nt(16)= float(totpt)
c$$$               nt(17)= float(ntrks)
c$$$               nt(18)= sngl(y0(KZ,IT))
c$$$               call hfn(4,nt)
c$$$            endif
c$$$#endif
            CHISQ = CHISQ + CH1 + 2.0D0*CH12 + CH2
+REP,TC_VERT,TCVRHX,1042.
      IF((ITER.GT.1).AND.(RCHISQ .GT. VHMCH)) THEN
+REP,TC_VERT,TCVRHX,1070-1073.
*     Check to see if it was an acceptable fit

      IF(RCHISQ .GT. VHGCH) IERR1 = 300
C      IF( PROB(SNGL(CHISQ),NDF).LT.VHCFL ) IERR1 = 300
+DEL,TC_VERT,TCVRHX,1081.
+REP,TC_VERT,TCVRHX,1091-1101.
*     Scale the (d(chi^2)/da^2)^-1 matrix by two, as given in the PDG
*     RPP, 1994, p. 1277, eqn. 17.26b
*     Covariance matrix = 2 * (d(chi^2)/da^2 matrix)^-1
*     At the same time, we are throwing away all other off diagonal
*     elements. This should be relatively safe, since the off diagonal
*     elements only tend to make the error ellipses have smaller area.

*     check values are within single precision bounds

      DO I=1,NPAR
         DO J=I,NPAR
            IF(ATGBA(I,J).LE.-1D20.OR.ATGBA(I,J).GE.1D20) THEN
               IERR1 = 700
               CALL ERRLOG(IERRL(6),
     &              '<TCVRHX> Matrix Element out of range (see log)')
               WRITE(LLOG,*) '<TCVRHX> ATGBA vertex error matrix:'
               DO II=1,NPAR
                  WRITE(LLOG,*) (ATGBA(II,JJ),JJ=1,NPAR)
               ENDDO
               GOTO 5000
            ENDIF
         ENDDO
      ENDDO

      Q(ITCVT+8)  = SNGL(2.D0*ATGBA(1,1))
      Q(ITCVT+9)  = SNGL(2.D0*ATGBA(2,1))
      Q(ITCVT+10) = SNGL(2.D0*ATGBA(2,2))
      Q(ITCVT+11) = SNGL(2.D0*ATGBA(3,1))
      Q(ITCVT+12) = SNGL(2.D0*ATGBA(3,2))
      Q(ITCVT+13) = SNGL(2.D0*ATGBA(3,3))

      Q(ITCVT+14) = SNGL(MIN(CHISQ,1D20))
+REP,TC_VERT,TCVRHX,1122-1127.
      MTCVP(4) = LENVP * NTRKS
*     ITCVT = LQ(LTCVX-IVRTX)
      IF (ITCVT1.LT.0) THEN
         ITCVT = XXTCVT(-ITCVT1)
      ELSE
         ITCVT = LQ(LTCVX-ITCVT1)
      ENDIF
      CALL MZLIFT(IXTCVX,ITCVP,ITCVT,-1,MTCVP,0)
* extra steps necessary in case of garbage collection
      IF (ITCVT1.LT.0) THEN
         ITCVT = XXTCVT(-ITCVT1)
      ELSE
         ITCVT = LQ(LTCVX-ITCVT1)
      ENDIF
      ITCVP = LQ(ITCVT-1)
+ADD,TC_VERT,TCVRHX,1134.
* extra steps necessary in case of garbage collection
      IF (ITCVT1.LT.0) THEN
         ITCVT = XXTCVT(-ITCVT1)
      ELSE
         ITCVT = LQ(LTCVX-ITCVT1)
      ENDIF
      ITCVP = LQ(ITCVT-1)
      ITCVH = LQ(ITCVP)
+DEL,TC_VERT,TCVRHX,1140-1222.
+REP,TC_VERT,TCVRHX,1238-1247.
*     Fill the new TCVH bank with the newly fitted helix parameters
*     converting the 3+3n parameters to 5n style

         IF (NTRKS.GT.1) THEN
            SGN1      = SGN(IT)
            ALPHA     = MIN(MAX(X(ALPHAX+IOFF),-1D6),1D6)
            TANLA     = MIN(MAX(X(TANLAX+IOFF),-1D2),1D2)
            INVALP    = 1.D0/ALPHA
+REP,TC_VERT,TCVRHX,1260-1265.
         ENDIF

         Q(ITCVH+1) = SNGL(MIN(MAX(R0,-100.D0),100.D0))

         Q(ITCVH+2) = SNGL(MIN(MAX(Z0,-50D0),50D0))
+REP,TC_VERT,TCVRHX,1290-1298.
*     Scale the (d(chi^2)/da^2)^-1 matrix by two, as given in the PDG
*     RPP, 1994, p. 1277, eqn. 17.26b.

         DO 6010 I=1,3
            DO 6020 J=1,3
               ATGBA(IOFF+I,IOFF+J) = 2.D0 * ATGBA(IOFF+I,IOFF+J)
 6020       CONTINUE
 6010    CONTINUE

*     Multiple Scattering Section -------------------

*     this code stolen from tcmsct(), and then heavily rewritten
*     Compute the sec^2 of lambda, (sec^2 = 1+tan^2).
*
         SECL2 = 1.0D0 + TANLA**2
         SECL  = SQRT(SECL2)
*
*           Compute the distance the particle has travelled in the JDC.
*             xl = chrg * (z_n - z_1) * [sec(lambda) / tan(lambda)]
*
         NPTSX = IQ(ITCTR+1)
         LYR0 = IQ(ITCTR+4)
         LYRN = IQ(ITCTR+5)
*
         IF((LYR0.LE.0).OR.(LYR0.GT.23)) LYR0 = 1
         IF((LYRN.LE.0).OR.(LYRN.GT.23)) LYRN = 23

*     The measurement error for N hits of average resolution
*     sigma_rf is
*
*                              [sigma_rf**2 * 720 ]
*     [sigma(alp)/alp]**2 = ---------------------------
*                           [ alpha * L**2 ]**2 * [N+4]
*
         SAA =  SIGRF / ( FLOAT(NPTSX+4)
     &        * ((RLYRTC(LYRN)-RLYRTC(LYR0))**2)**2)
*
         IF (SAA.GT.ATGBA(IOFF+ALPHAX,IOFF+ALPHAX)) THEN
            ATGBA(IOFF+ALPHAX,IOFF+ALPHAX)= DBLE(SAA)
         ENDIF
*
*     The measurement error of tan(lambda) for N hits, and
*     track length L is given as:
*
*                                          sig_z**2   12(N-1)
*    [sigma_tan(lam)]**2 = sec**2(lambda)* -------- * -------
*                                          del_r**2   N(N+1)
*
         SAA = SIGZ*(1.D0+TANLA**2)*FLOAT(NPTSX-1)
     &        /      (FLOAT(NPTSX*(NPTSX+1))
     &        *      (RLYRTC(LYRN)-RLYRTC(LYR0))**2)
*
         IF (SAA.GT.ATGBA(IOFF+TANLAX,IOFF+TANLAX)) THEN
            ATGBA(IOFF+TANLAX,IOFF+TANLAX) = DBLE(SAA)
         ENDIF
*
*     XLEN is the path length through the JDC
*
         XLEN = SECL * (RLYRTC(LYRN) - RLYRTC(LYR0))
         IF(XLEN .LE. 0) XLEN = 1.0
*
*     Compute the term [alpha/eBsec(lambda)]**2, (one over the
*     total momentum squared), then compute the total momentum
*     and the total energy squared.
*
         PTOT = (PP*SECL)**2
         IF (PTOT.LT.1E-4) PTOT=1E-4
         PINV = 1.D0/PTOT
         ETOT = PTOT + SQMPI
*
*     Multiple scattering eqn, Review of Part Prop.1994  eqn #10.8
*
*         13.6 MeV/c
*     t = ------------ * Z * sqrt(x/x0) * (1 + 0.038 ln (x/x0))
*          p * beta
*
*     The units of the 13.6 are MeV/c
*
         ALM  = (13.6**2)* ETOT * PINV * PINV
*
*     Now calculate the multiple-scattering contribution to the
*     covariance matrix.
*
*     These are GAS scattering corrections. The hard stuff corrections
*     are easier.
*     Gluckstern, p 384, eq (21)
*
*     MS errors for circle, thus rho/phi projection
*     <c^2>    =  K Cn / L
*     <c theta> = -K Dn
*     <theta^2> = K L En
*
*     MS errors for straight line, thus rho/zed projection
*     <theta(SL)^2> = 1/3 K L
*
*     Translation:
*     c         = curvature (-alpha*sgn)
*     theta     = azimuthal psi/phi
*     K         = MS angle per distance of scatterer
*     L         = XLEN
*     theta(SL) = lambda
*
*     d(theta(SL))^2 * sec^4 = d(tanlambda)^2
*
*     Note that K, Cn, Dn and En are independent of orientation,
*     inside the JDC gas.  Thus the only thing that should depend on
*     track angle (tan lamda) should be L = XLEN, which already has
*     the cosecant correction. Therefore, there shouldn't be all
*     these extra factors of tan lambxda.

*     GASXT is (gas radiation length)^-1

         ALMGAS = ALM * GASXT * (1+0.038*LOG(GASXT*XLEN))
         S33 =  ALMGAS * CN / XLEN
         S44 =  ALMGAS * 0.3333 * XLEN * SECL2*SECL2
         S55 =  ALMGAS * EN * XLEN
         S35 =  ALMGAS * DN * SGN1

*     These are pipe+pwc+jdc inner wall corrections, only
*     needed if vertex is inside of the inner jdc wall.
*     Note that we are lumping all the corrections together
*     even if the radius is somewhere in between.
*     Please find the correct vacuum pipe radius, this is a guess

*     Include the errors of the radius to the radius

         RR = SQRT(MAX(X(1)**2+X(2)**2-ATGBA(1,1)-ATGBA(2,2),0.0D0))

*     We have 6 significant sources of matter. We'll see how many layers
*     the track needs to pass through, if we know where the vertex is.

*     1 = hydrogen
*     2 = mylar
*     3 = Al vacuum pipe
*     4 = PWC1
*     5 = PWC2
*     6 = JDC inner wall

         DO 6100 IRHO=1,6
            IF (RR.LE.SCTRAD(IRHO)) THEN
*     radial error
               S11 = ALM*SECL*SECL2*R2X(IRHO)
*     zed error
               S22 = ALM*SECL*SECL2*R2X(IRHO)
*     tan lambda error
               S44 = S44 + ALM*SECL2*SECL2*SECL*R0X(IRHO)
*     psi/phi error
               S55 = S55 + ALM*SECL*R0X(IRHO)
* break out of loop
               GOTO 6200
            ENDIF
 6100    CONTINUE
 6200    CONTINUE

*     finally, add in the Multiple Scattering errors to the helix
*     parameters.
         ATGBA(IOFF+1,IOFF+1) = ATGBA(IOFF+1,IOFF+1) + S33
         ATGBA(IOFF+2,IOFF+2) = ATGBA(IOFF+2,IOFF+2) + S44
         ATGBA(IOFF+3,IOFF+3) = ATGBA(IOFF+3,IOFF+3) + S55
         ATGBA(IOFF+1,IOFF+3) = ATGBA(IOFF+1,IOFF+3) + S35
         ATGBA(IOFF+3,IOFF+1) = ATGBA(IOFF+3,IOFF+1) + S35

*     Also add the Multiple Scattering errors to the vertex position

         Q(ITCVT+8)  = Q(ITCVT+8)  + S11
         Q(ITCVT+10) = Q(ITCVT+10) + S11
         Q(ITCVT+13) = Q(ITCVT+13) + S22

*     this is the transformation matrix from (alpha,tanla,phi) to (p1,p2,p3)
*
*     dpx/dalpha dpy/dalpha dpz/dalpha
*     dpx/dtanla dpy/dtanla dpz/dtanla
*     dpx/dphi   dpy/dphi   dpz/dphi
+REP,TC_VERT,TCVRHX,1304-1312.
*     row 2, d tan lam
         DPDA(2,1) = 0.D0
         DPDA(2,2) = 0.D0
         DPDA(2,3) = PP
*     row 3, d phi
         DPDA(3,1) = PP*COSB*SGN1
         DPDA(3,2) = PP*SINB*SGN1
         DPDA(3,3) = 0.D0

* matrix arithmetic (3x3):
*              dp T           dp
* mom error = (--)  . atgba . --
*              da             da
*
+REP,TC_VERT,TCVRHX,1336-1341.
*
***SPP(I,J) could be too much high which causes the code to crash...a limit is set then...
*
         DO I=1,3
            DO J=1,3
               IF(DABS(SPP(I,J)).GT.1.D+10) SPP(I,J)=1.D+10
            ENDDO
         ENDDO
*
*     | 8        |
*     | 9  10    |
*     | 11 12 13 |
*     protect against large values!

         Q(ITCVP+ 8) = SNGL(MIN(SPP(1,1),1.D8))
         Q(ITCVP+ 9) = SNGL(MAX(MIN(SPP(1,2),1.D8),-1.D8))
         Q(ITCVP+10) = SNGL(MIN(SPP(2,2),1.D8))
         Q(ITCVP+11) = SNGL(MAX(MIN(SPP(1,3),1.D8),-1.D8))
         Q(ITCVP+12) = SNGL(MAX(MIN(SPP(2,3),1.D8),-1.D8))
         Q(ITCVP+13) = SNGL(MIN(SPP(3,3),1.D8))


+REP,TC_VERT,TCVRHX,1369-1370.
* BUG FIX, I did it this way for readability, but NTRKS is readonly
CMTL         NTRKS = 0
CMTL         IQ(ITCVT+1) = NTRKS
         IQ(ITCVT+1) = 0
+DEL,TC_VERT,TCVRHX,1376.
+ADD,TC_VERT,TCVRHZ,0.
*CMZU:  2.01/01 13/03/96  14.43.37  by  F.-H.Heinsius
+REP,TC_VERT,TCVRHZ,63-64.
                  IF (IQ(ITCVT+4).GT.0) THEN
                     CHI = Q(ITCVT+14)/FLOAT(IQ(ITCVT+4))
+ADD,TC_VERT,TCVPCI,*.
&DECK,TCVPCI.
*CMZU:  2.01/05 02/04/96  15.41.25  by  F.-H.Heinsius
*-- Author :
      SUBROUTINE TCVPCI(NAME,VALUE,DEFAULT)
&SEQ,IMPNONE.
&SEQ,CBUNIT.
      CHARACTER*(*) NAME
      INTEGER VALUE,DEFAULT
      IF ( VALUE.NE.DEFAULT) THEN
         WRITE(LLOG,295) NAME,VALUE,DEFAULT
      ELSE
         WRITE(LLOG,294) NAME,VALUE
      ENDIF
 294  FORMAT(' **TCVRHX: ',A,' is default ',I10)
 295  FORMAT(' **TCVRHX: ',A,' was ',I10,', is now ',I10)
 
      END
+ADD,TC_VERT,TCVPCR,*.
&DECK,TCVPCR.
*CMZU:  2.01/13 14/03/97  15.40.29  by  M.T. Lakata
*CMZU:  2.01/05 02/04/96  15.41.25  by  F.-H.Heinsius
*-- Author :
      SUBROUTINE TCVPCR(NAME,VALUE,DEFAULT)
&SEQ,IMPNONE.
&SEQ,CBUNIT.
      CHARACTER*(*) NAME
      REAL VALUE,DEFAULT
      IF ( VALUE.NE.DEFAULT) THEN
         WRITE(LLOG,295) NAME,DEFAULT,VALUE
      ELSE
         WRITE(LLOG,294) NAME,VALUE
      ENDIF
 294  FORMAT(' **TCVRHX: ',A,' is default ',F12.5)
 295  FORMAT(' **TCVRHX: ',A,' was ',F12.5,', is now ',F12.5)
 
      END
+REP,TC_SVTX,RSRVTX,0.
&DECK,RSRVTX.
*CMZU:  2.01/12 03/03/97  15.23.15  by  Rafik Ouared
+REP,TC_SVTX,TSTVTX,0.
&DECK,TSTVTX.
*CMZU:  2.01/12 03/03/97  15.23.42  by  Rafik Ouared
+ADD,TC_SVTX,TSVPOS,0.
*CMZU:  2.01/02 14/03/96  13.27.47  by  F.-H.Heinsius
+ADD,TC_SVTX,TSVPOS,20.
*        Modified   :   14,March,1996 R.Ouared
*                       MTPCH replaced by MTSV1 for both BP amplitudes
*                       and 'hits' banks and by MTSV2 for 'clusters'
*                       bank. The number of matched clusters set to 0
*                       in tvxt bank (IQ(LTVXT+6)=0).
+REP,TC_SVTX,TSVPOS,38-42.
*        K     INDEX FOR CLUSTER NUMBER
*        ND    Number of data words in the RPWC bank.
*        NWIRE  The number of wires in CLUSTER.
+REP,TC_SVTX,TSVPOS,49.
      INTEGER  I,II,ND,NWIRE,IWIR1,J,JWIR1,KK,NCONT,L
+REP,TC_SVTX,TSVPOS,59-62.
*        Pointer for the processed data.
*        ITVTX
*
      INTEGER  ITVTX
+REP,TC_SVTX,TSVPOS,171-172.
      MTSV1(4) = 16 * LENBP
      CALL MZLIFT(IXTVXT,ITVTX,LTVXT,-3,MTSV1,0)
+REP,TC_SVTX,TSVPOS,323-324.
      MTSV1(4) = II * LENSV
      CALL MZLIFT(IXTVXT,ITVTX,LTVXT,-1,MTSV1,0)
+ADD,TC_SVTX,TSVPOS,364.
      IQ(LTVXT+6) = 0                          !# of matched clusters.
+REP,TC_SVTX,TSVPOS,389-390.
      MTSV2(4) = II * LENCLS
      CALL MZLIFT(IXTVXT,ITVTX,LTVXT,-2,MTSV2,0)
+ADD,TC_UTIL,TCPUSH,0.
*CMZU:  2.01/05 02/04/96  19.14.40  by  F.-H.Heinsius
+REP,TC_UTIL,TCPUSH,10.
      INTEGER VALUE,COUNT,IDATA
+ADD,TC_UTIL,TCPOP,0.
*CMZU:  2.01/05 02/04/96  19.13.56  by  F.-H.Heinsius
+REP,TC_UTIL,TCPOP,10.
      INTEGER VALUE,COUNT,IDATA