banaspi.f

C----------------------------------------------------------------------
C
C  This file contains:
C
C   Subroutine BANASPI (ITRACK,XMAS,IRC,ICHG,BCHISQ,XI0,TKP,ERR,IVERSION)
C
C
C     BANASPI is used to recover the uncorrected
C     error matrix for a charged track (ICECREAM) and apply
C     corrections to the track parameters and error matrix
C     corresponding to a particular mass hypothesis (FUDGE).
C     A BMFIT is then performed.
C     DUE to the lack of trklist message,
C     the uncorrected error matrix for a charged track can not
c     be recovered for iversion less than 8.we just keep this
c     subroutine for doing bmfit with the error matrix from
c     trklist for iversion less than 8.For iversion.ge.8,one alwayes
c     use banaspig.           H.B.Lan,Jan. 14,1994
C
C     Input:
C        I*4  ITRACK      TRKLST track number
C        R*4  XMAS        desired mass hypothesis
C        I    IVERSION    THE RECONSTRUCTED VERSION     ! LANHB
C     Output:
C        I*4  IRC         Return code:
C                          0  Normal completion
C                          1  Track is not charged!
C                          2  See discussion of BANASPIG below.
C                          3  Track does not have MFIT=2.
C                          4  BMFIT failure.
C
C        I*4  ICHG        Track charge
C        R*4  BCHISQ      chi-squared from BMFIT
C        R*4  XI0         xi after FUDGE (before BMFIT)
C                           = SIGNED distance of closest
C                               approach to BEAM AXIS (meters)
C        R*4  TKP(3)      phi, 1/Pxy, tan(lambda)
C        R*4  ERR(6)      1st 6 elements of error matrix
C
C       If the geometric (undamaged) error matrix is available
C       in your data set, you may want to use the entry point
C       BANASPIG.  This is faster and gives more accurate results
C       since the error matrix does not need to undergo the
C       approximate ICECREAM process.  If you call BANASPIG
C       and the geometric error matrix is not available, you
C       will get IRC=2.  BANASPIG has the same parameters
C       as BANASPI.
C
C                                    30 Sep 87  SRW
C
C
C   Make psdobmft also available for iversion less than 8.
C                                    Jan. 12,1994
c
c       A seperate subroutine at BES,which are useful only for
c       iversion less than 8
C                                     Jan. 14,1994
C----------------------------------------------------------------------

banaspig.f

C----------------------------------------------------------------------
C
C  This file contains:
C
C   Subroutine BANASPIG (ITRACK,XMAS,IRC,ICHG,BCHISQ,XI0,TKP,ERR,IVERSION)
C
C       If the geometric (undamaged) error matrix is available
C       in your data set, you may want to use the entry point
C       BANASPIG.  This is faster and gives more accurate results
C       since the error matrix does not need to undergo the
C       approximate ICECREAM process.  If you call BANASPIG
C       and the geometric error matrix is not available, you
C       will get IRC=2. for iversion.ge.8,the geometric matrix is
c       available.             H.B.Lan  Jan. 14,1994
C     Input:
C        I*4  ITRACK      TRKLST track number
C        R*4  XMAS        desired mass hypothesis
C        I    IVERSION    THE RECONSTRUCTED VERSION     ! LANHB
C     Output:
C        I*4  IRC         Return code:
C                          0  Normal completion
C                          1  Track is not charged!
C                          2  See discussion of BANASPIG below.
C                          3  Track does not have MFIT=2.
C                          4  BMFIT failure.
C
C        I*4  ICHG        Track charge
C        R*4  BCHISQ      chi-squared from BMFIT
C        R*4  XI0         xi after FUDGE (before BMFIT)
C                           = SIGNED distance of closest
C                               approach to BEAM AXIS (meters)
C        R*4  TKP(3)      phi, 1/Pxy, tan(lambda)
C        R*4  ERR(6)      1st 6 elements of error matrix
C
C-----------------------------------------------------------------------

bmfit.f

C-----------------------------------------------------------------------
C
C  This file contains:
C
C   Subroutine BMFIT(ITRACK,B,BCHISQ)
C
C
C   NOTE: This routine is provided to simulate the operation of
C  the original Bmfit routine.  This routine uses the Trklst
C  Subtype 3 track parameters and error matrix as input to the
C  beam fit.  It is preferable to first make corrections for
C  dE/dx and multiple scattering, and then to perform the beam fit.
C  See BANASPI (in Bsplit Fortran).
C
C   Bmfit (Beam fit) is intended to replace the TRKLST Subtype 8
C  Vertex Constrained fit.  The Beam fit effectively adds the
C  Interaction Region as a measured track point.  (The fit is NOT
C  really made at the DC hit level, however.)  The measured IR position
C  and size are taken from TRKLST Vertex Subtype 3 (=BEAM CROSSx).
C
C  BMFIT uses the NAPL subroutine DPPFCS, therefore one needs to have the
C  NAPL TXTLIB accessed.
C
C  BMFIT calls SWIMIT to move the track parameters and their errors along
C  the track to the point of closest approach in the x-y plane to the beam
C  position.
C
C  Then the errors in the beam position are used in performing a least
C  squares calculation to minimize this distance of the track to the beam.
C
C  New track parameters and their error matrix are returned.
C
C  Input:
C    I*4  ITRACK   TRKLST track number
C
C  Output:
C    R*8  B(5,6)   Results of fit:
C           New track parameters given by B(1-5,1):
C                  B(1,1) = phi
C                  B(2,1) = kappa
C                  B(3,1) = tan(lambda)
C                  B(4,1) = xi
C                  B(5,1) = eta
C           New variance matrix given by B(1-5,2-6);
C           e.g.:  B(1,2) = variance of phi
C                  B(3,6) or B(5,4) = covariance of tan(lambda) and eta
C                  etc.
C
C    R*4  BCHISQ   Chi-squared of fit;
C               should follow distribution for one degree of freedom;
C               set to -1 if unable to perform fit.
C
C
C    Original Version ....................... Larry Parrish  19-Feb-1987
C
C------------------------------------------------------------------------
C
C    Placed the innards in SUBROUTINE BMFIT0 to permit use of the
C    fitter on dE/dx-corrected tracks.
C                     ................................. SRW  24-Sep-1987
C
C-----------------------------------------------------------------------

bmfit0.f

C-----------------------------------------------------------------------
C
C  This file contains:
C
C   Subroutine BMFIT0(ICHG,TKPAR,ERROR,B,BCHISQ,XI0)
C
C
C   Bmfit0 is the heart of L. Parrish's original Bmfit.
C  For maximum rigor, obtain the corrected track parameters
C  and error matrix (e.g. using other Telesis routines) and
C  call this routine.  Subroutine BANASPI is provided to
C  perform the full procedure given the TRKLST track number
C  and mass hypothesis.
C
C  See Bmfit Mortran for more information about the Beam fit.
C
C  Input:
C    I*4  ICHG        track charge
C    R*4  TKPAR(6)    phi,k,s,xm,ym,eta
C    R*4  ERROR(15)   error matrix
C
C  Output:
C    R*4  TKPAR(6)    Fitted phi,k,s,xm,ym,eta
C    R*4  ERROR(15)   New error matrix
C    R*8  B           5 x 6 double precision real matrix,
C                       contains results of fit;
C           New track parameters given by B(1-5,1):
C                  B(1,1) = phi
C                  B(2,1) = kappa
C                  B(3,1) = tan(lambda)
C                  B(4,1) = xi
C                  B(5,1) = eta
C           New variance matrix given by B(1-5,2-6);
C           e.g.:  B(1,2) = variance of phi
C                  B(3,6) or B(5,4) = covariance of tan(lambda) and eta
C                  etc.
C
C    R*4  BCHISQ      chi-squared of fit (1 degree of freedom);
C                      = -1 if unable to perform fit, or
C                           if track charge is flipped by fit.
C
C    R*4  XI0         xi [= (xm-xir)*sin(phi) - (ym-yir)*cos(phi)
C                         = signed impact distance] measured from
C                         the BEAM AXIS, BEFORE THE FIT.
C                         [Notice that xi0 can be less than zero.
C                          Cut on abs(xi0) instead of RM.]
C
C                                   SRW  21 Apr 1988
C
C-----------------------------------------------------------------------

dmpsqdl.f

C-----------------------------------------------------------------------
C
C  This file contains:
C
C   Subroutine DMPSQDL(N_PART,MASS,P_SPEC,N_RES,M_RES,P_MEM,MEM,RC,CH 
C                      ISQ,P_FIT)
C
C
C   DMPSQDL gives a printed dump of the input, output and values of
C the Squid Common block.
C                                        John Brown  85/09/02
C
C-----------------------------------------------------------------------

fudge.f

C-----------------------------------------------------------------------
C
C  This file contains:
C
C   SUBROUTINE     FUDGE(XMAS,TKPAR,ICHG,BC,PM,SECLM,IVERSION)
C
C   The SUBROUTINE Fudge is provided to make dE/dx corrections
C  to a track's momentum when the error matrix is not of interest.
C  (Of course this is significantly faster than SUPERFUG.)
C  The input parameters are XMAS, TKPAR, ICHG, BC, PM, and SECLM,
C  as defined above.  The corrected track parameters are output
C  in TKPAR.
C  Input:
C     R*4   XMAS        Mass hypothesis (GeV/c**2);
C     R*4   TKPAR(6)    Vector of track parameters:
C                         phi, k (=1/Pxy), s (=tanlam), xm, ym, zm;
C     I*4   ICHG        Charge;
C     R*4   BC          = 0.029979 * B field at (xm,ym,zm);
C     R*4   PM(5)       Vector containing Px, Py, Pz, Pxyz, Pxy;
C     R*4   SECLM(4)    Ith element = sec(lambda)**I;
C     R*4   ERROR(15)   Error matrix (taken from TRK(I3+20) et seq.);
C     R*4   RM          = sqrt(xm**2 + ym**2);
C     I*4   MFIT        = Least-significant byte of MFITWD;
C     I*4   MZFIT       = Next-to-least-significant byte of MFITWD;
C     I*4   IFIT(2)     R-phi layer maps = ITRK(I3+35), ITRK(I3+36);
C     I*4   CD          CD CODE = ITRK(I3+40).
C
C Output:
C           TKPAR(6)    Corrected track parameters;
C           ERROR(15)   Corrected error matrix.
C
C   The following track parameters are corrected: phi, k, xm, and ym.
C
C
C
C--------------------------------------------------------------------C
C
C  Miniature set of COULOM-like constants:
C
C  U        = exponent in D. Hitlin formula
C  XMCS     = (0.0141 GeV)**2 for Multiple Scattering
C  NSCATR   = number of scatterers under consideration
C
C               In the following arrays, the scatterers
C               are numbered as follows:
C                  1. Beampipe
C                  2. Little Chamber
C                  3. Inner wall of Main DC
C                  4. Gas and wires of DC layers 2-5
C
C  RSCATR   = radius of mid-scatterer from beam axis
C  DPMIN    = d(momentum) for each scatterer (from Hitlin)
C               for a minimum-ionizing particle
C  ALPHAT   = total alpha*t for each scatterer, determined
C               from COULOM COMMON.  These numbers are
C               not necessarily consistent with the DPMIN
C               values, but it's not worth worrying about.
C  RADL     = number of radiation lengths in each scatterer
C
C
C                   A SEPERATE SUBROUTINE AT BES
C                             H.B.LAN   JAN.14,1994
C--------------------------------------------------------------------C

ggfit.f

C-----------------------------------------------------------------------
C
C  This file contains:
C
C   Subroutine GGFIT(ITK1,ITK2,GGMASS,IRC,RAWMGG,RAWCOS,  CHISQ,FITCOS
C                    *,FITPG1,FITPG2,FITPGG)
C
C
C   Ggfit is used to perform a 1C kinematic fit to pi0 -> gamma gamma
C  or eta -> gamma gamma (or anything -> gamma gamma).  Of course
C  one can always use Squidjr directly, but this routine does the
C  fit and a few other common calculations with a single call.
C  The user specifies the Trklst track numbers of the two photon
C  candidates, and the desired gamma-gamma mass constraint.  Out
C  comes the unfit mass and decay angle, and the fit chi-squared,
C  decay angle, and four-vectors.
C
C   Needless to say, one must check the Ggfit return code before
C  using the results.
C
C   The `decay angle' is defined to be the angle in the gamma-gamma
C  center-of-mass system between the laboratory frame and one of
C  the photons.  [Since we take abs(cos(decay angle)) it doesn't
C  matter WHICH photon is used.]  In fact
C
C           abs(cos(decay angle)) = abs(E1 - E2)/Pgg
C
C  where E1 and E2 are the lab gamma energies and Pgg is the lab
C  momentum of the gamma-gamma system.  The proof is left as an
C  exercise.
C
C  Input:
C     I*4   ITK1,ITK2    Trklst track numbers
C     R*4   GGMASS       Mass (GeV) for gamma-gamma mass constraint
C
C  Output:
C     I*4   IRC          Return code
C                         = 0  Normal completion;
C                           1  Specified track(s) absent or not neutral;
C                           2  Squidjr fit failed.
C     R*4   RAWMGG       Raw gamma-gamma mass (GeV)
C     R*4   RAWCOS       abs(cos(decay angle)) using raw gammas
C     R*4   CHISQ        Fit chi-squared (one degree of freedom)
C     R*4   FITCOS       abs(cos(decay angle)) using fit gammas
C     R*4   FITPG1(4)    Four-momentum of gamma 1 after fit
C     R*4   FITPG2(4)    Four-momentum of gamma 2 after fit
C     R*4   FITPGG(4)    Four-momentum of gamma gamma after fit
C
C
C                            Steve Wasserbaech  15 April 1988
C
C
C-----------------------------------------------------------------------

kfinit.f

C-----------------------------------------------------------------------
C   Subroutine KFINIT (NPART,MASS,PSPEC,NRES,MRES,PMEM,MEM,SQCHI,PFIT)
C
C Set up parameters for constrained fit:
C    Load X  from TRKLST, set ETOT, NC, NVAR. 
C
C    CASE=1  : No missing particles or components
C        =2  : Missing 3-momenta, missing mass known
C        =3  : Missing energy for neutral particle
C
C-----------------------------------------------------------------------

kfinitjr.f

C-----------------------------------------------------------------------
C
C  This file contains:
C
C   Subroutine KFINITJR(N_PART,MASS,P_SPEC,N_RES,M_RES,P_MEM,MEM,RC,
C                      *CHISQ,P_FIT)
C
C
C   Kfinitjr (an internal Squidjr routine) makes a few preparations
C before starting on the fit.
C
C   In the old Squid this routine was a bit more complicated.  Most
C of its responsibilities have been moved to the pre-squid routines.
C
C   The parameters of this routine are the same as those of
C Squidjr itself.
C
C SRW 19 May 88: Add initialization of chisq in case fit is satisfied
C on zeroth iteration.
C    A seperate subroutine at BES     H.B.Lan ,Jan. 14 ,1994
C-----------------------------------------------------------------------

kfinitl.f

C-----------------------------------------------------------------------
C
C  This file contains:
C
C   Subroutine KFINITL(N_PART,MASS,P_SPEC,N_RES,M_RES,P_MEM,
C                     *MEM,RC,CHISQ,P_FIT)
C
C
C   KFINITL (an internal Newsquid routine) makes a few preparations
C before starting on the fit.THIS SUBROUTINE JSUT INITIALIZE SOME
C QUANTITIES.
C   THE REAL CALCULATION OF DIFFERENT VARIATIES WAS CARRIED OUT IN
C THE SUBROUTINE KFUPDTL             H.B.LAN,    JAN. 7,1994
C
C   In the old Squid this routine was a bit more complicated.  Most
C of its responsibilities have been moved to the pre-squid routines.
C
C   The parameters of this routine are the same as those of
C Newsquid itself.
C
C SRW 19 May 88: Add initialization of chisq in case fit is satisfied
C on zeroth iteration.
C
C       MOVE OUT THE ENTRY POINT KFUPDTL TO A SEPERATE SUBROUTINE,
C       IN ORDER TO FIX THE MACHINE-RELATED BUG FOUND BY B.LOWERY
C                                             H.B.LAN
C                                             JAN. 7,1994
C
C
C-----------------------------------------------------------------------

kfupdt.f

C-----------------------------------------------------------------------
C        SUBROUTINE KFUPDT(NPART,MASS,PSPEC,NRES,MRES
C                          ,PMEM,MEM,SQCHI,PFIT)
C
C calculate constraints and derivatives from X
C-----------------------------------------------------------------------

kfupdtjr.f

C-----------------------------------------------------------------------
C
C  This file contains:
C
C   Subroutine KFUPDTJR(N_PART,MASS,P_SPEC,N_RES,M_RES,P_MEM,MEM,RC,
C                      *CHISQ,P_FIT)
C
C   Kfupdtjr  calculate constraints and the derivatives from X
C   In the old Squid this routine was a bit more complicated.  Most
C of its responsibilities have been moved to the pre-squid routines.
C
C   The parameters of this routine are the same as those of
C Squidjr itself.
C
C SRW 19 May 88: Add initialization of chisq in case fit is satisfied
C on zeroth iteration.
C             A seperate subroutine at BES
c
c                                           H.B.Lan, Jan.14,1994
C-----------------------------------------------------------------------

kfupdtl.f

C-----------------------------------------------------------------------
C
C  This file contains:
C
C   Subroutine KFUPDTL(N_PART,MASS,P_SPEC,N_RES,M_RES,P_MEM,
C                     *MEM,RC,CHISQ,P_FIT)
C
C
C
C   In the old Squid this routine was a bit more complicated.  Most
C of its responsibilities have been moved to the pre-squid routines.
C
C   The parameters of this routine are the same as those of
C Newsquid itself.
C
C SRW 19 May 88: Add initialization of chisq in case fit is satisfied
C on zeroth iteration.
C
C----------------------------------------------------------------------
C  Calculate constraints and derivatives from X.
C
C-----------------------------------------------------------------------
C
C   IT IS EXTRACTED FROM ENTRY POINT IN SUBROUTINE
C   KFINITL ,IN ORDER TO FIX THE MACHINE-
C   RELATED BUG FOUND BY B.LOWERY
C
C                                    H.B.LAN
C                                    JAN. 7,1994
C   fix one bug related to the case of photon missing energy
C                                    H.B.Lan
C                                    14,July,1994
C-----------------------------------------------------------------------

klamcopy.f

C-----------------------------------------------------------------------
C       SUBROUTINE KLAMCOPY(ITK1,ITK2,Q1,Q2,MASS1,MASS2,TK1,TK2,
C     * ER1,ER2,IRC,TKS1,TKS2,ERS1,ERS2,CRXY,ZSWUM,P4PAR,MPAR,
C     * PXYZ,CHI2XY,CHI2Z,LXY,CT,CHISQFIT,FLAGFIT,iversion)
C 
C Called by KlamCORR
C-----------------------------------------------------------------------


klamcorr.f

C-----------------------------------------------------------------------
C
C   This file contains:
C
C    Subroutine KlamCORR(ITK1,ITK2,MASS1,MASS2,IRC,TKS1,TKS2,ERS1,ERS2,
C    *CRXY,ZSWUM,P4PAR,MPAR,PXYZ,CHI2XY,CHI2Z,LXY,CT,CHISQFIT,
C    * FLAGFIT,IVERSION)
C
C
C    KLAMCORR does the corrections for
C  dE/dx and multiple scattering are actually computed.  Use this
C  routine when the dE/dx-corrections blocks are not available in
C  Trklst track subtype 3 for the pion and/or proton hypothesis.
C
C   The heart of the calculation is placed in Subroutine KlamCOPY,
C  allowing the user to save cpu time if the same charged tracks
C  are to be used in several combinations.  The user calls
C  SUPERSD or SUPERSDg for each track to perform the
C  dE/dx corrections, then loops over combinations calling
C  KlamCOPY.
C
C  Input:
C     I*4   ITK1,ITK2    Trklst track numbers
C     R*4   MASS1,MASS2  Masses (GeV) of tracks 1 and 2
C     I     IVERSION     RECONSTRUCTED VERSION   !LANHB
C  Output:
C     I*4   IRC        Return code
C                       = 0  Normal completion;
CLHB                    NOT EQUAL 0 OTHERWISE
C                         1  Specified track(s) not charged;
C                         3  Specified track(s) not helix fit;
C                         4  Tracks are not oppositely charged;
C                         5  Mass1, Mass2 are not consistent
C                                with pi pi or p pi;
C                         6  Parent has Pxy = 0 before swimming;
C                         7  Daughter has Pxy = 0;
C                         8  Track circles do not intersect;
C                         9  Two good crossing points found;
C                        10  Zero good crossing points found;
C                        11  Parent has Pxy = 0 after swimming;
C                        12  Arithmetic error computing uncertainty
C                                in crossing point position;
C
C                         For IRC=0 all outputs are filled.
C
C                         Fatal return codes are 1,3,5,6,7,8,11,12;
C                       none of the other outputs are filled.
C
C                         For IRC=9 the closer crossing point to the
C                       interaction point is chosen and all outputs
C                       are filled.
C
C                         For IRC=4 or 10 no swimming is done.
C                       TKS1, TKS2, ERS1, ERS2, P4PAR, MPAR, PXYZ
C                       are filled using the dE/dx corrected (unswum)
C                       track parameters.  CRXY, ZSWUM, CHI2XY,
C                       CHI2Z, LXY, CT are not filled.
C
C     R*4   TKS1(3)    Track 1 swum parameters (phi,1/Pxy,tanlam)
C     R*4   TKS2(3)    Track 2 swum parameters
C     R*4   ERS1(6)    Track 1 swum error matrix (usual sequence)
C     R*4   ERS2(6)    Track 2 swum error matrix
C     R*4   CRXY(2)    (x,y) position of chosen crossing poing (meters)
C     R*4   ZSWUM(2)   z position of tracks 1 and 2 at crossing point (m)
C     R*4   P4PAR(4)   Parent (Px,Py,Pz,E) after swimming (GeV)
C     R*4   MPAR       Parent invariant mass (GeV)
C     R*4   PXYZ       Magnitude of parent three-momentum (GeV)
C     R*4   CHI2XY     Chi-squared (1 dof) of xy alignment of parent
C                       momentum with line from interaction point
C                       to crossing point.  See Mark III Memo
C                       referenced above.
C     R*4   CHI2Z      Chi-squared (1 dof) for the hypothesis that
C                       the two tracks have equal z coordinates at
C                       the crossing point.  **WARNING**: This
C                       quantity would have no meaning if one or
C                       both of the tracks has mzfit .ne. 3, so
C                       CHI2Z = -1.0 is returned in this case.
C     R*4   LXY        Best fit value of lab xy decay distance (m)
C     R*4   CT         c * proper decay time (m), derived from LXY,
C                       P4PAR, and the exact K0 or Lambda mass.
C
C
C                            Steve Wasserbaech  11 April 1988
CLHB---------FLAGFIT=.TRUE. DO FITVEE--AND CHISQFIT CORRESPONDES TO
CLHB---------THE CHI SQUARE------------------------------------------
C                               HUIBIN LAN 14 JAN 1993
C
C------------------------------------------------------------------------

klams.f

C-----------------------------------------------------------------------
C
C   This file contains:
C
C    Subroutine Klams(ITK1,ITK2,MASS1,MASS2,  IRC,TKS1,TKS2,ERS1,ERS2,
C     *CRXY,ZSWUM,  P4PAR,MPAR,PXYZ,CHI2XY,CHI2Z,LXY,CT,IVERSION)
C
C
C   KLAMS is used to select candidates for K0 -> pi+ pi- or
C  Lambda -> p pi-.  The user inputs the Trklst track numbers of two
C  charged tracks to be considered, and the assumed masses of these
C  tracks.  The routine looks up the dE/dx-corrected track parameters
C  in Trklst subtype 3 and determines the crossing points of the track
C  circles projected onto the xy-plane.  The correct crossing point
C  is chosen and the track parameters and error matrices are swum
C  to that point.  This information can easily be loaded
C  into Newsquid or Squidjr using Psexternal.  The pi pi or p pi
C  four-momentum and invariant mass is computed.  The decay
C  length is determined using the Decay Length Formalism
C  described in Mark III Memo 11/86-4 by Wasserbaech and Bolton.
C  Quantities are also returned which measure consistency of
C  the measured parent (K0 or Lambda) momentum direction with that
C  of the segment joining the interaction point and the crossing
C  point, and consistency of the z coordinates of the two
C  tracks at their crossing point (see below).
C
C   The heart of the calculation is placed in Subroutine Klams0.
C
C   Needless to say, one must check the Klams return code before
C  using the results.
C
C  Input:
C     I*4   ITK1,ITK2    Trklst track numbers
C     R*4   MASS1,MASS2  Masses (GeV) of tracks 1 and 2
C
C  Output:
C     I*4   IRC        Return code
C                       = 0  Normal completion;
C                         1  Specified track(s) not charged;
C                         2  Trklst dE/dx correction block not filled;
C                         4  Tracks are not oppositely charged;
C                         5  Mass1, Mass2 are not consistent
C                                with pi pi or p pi;
C                         6  Parent has Pxy = 0 before swimming;
C                         7  Daughter has Pxy = 0;
C                         8  Track circles do not intersect;
C                         9  Two good crossing points found;
C                        10  Zero good crossing points found;
C                        11  Parent has Pxy = 0 after swimming;
C                        12  Arithmetic error computing uncertainty
C                                in crossing point position;
C
C                         For IRC=0 all outputs are filled.
C
C                         Fatal return codes are 1,3,5,6,7,8,11,12;
C                       none of the other outputs are filled.
C
C                         For IRC=9 the closer crossing point to the
C                       interaction point is chosen and all outputs
C                       are filled.
C
C                         For IRC=4 or 10 no swimming is done.
C                       TKS1, TKS2, ERS1, ERS2, P4PAR, MPAR, PXYZ
C                       are filled using the dE/dx corrected (unswum)
C                       track parameters.  CRXY, ZSWUM, CHI2XY,
C                       CHI2Z, LXY, CT are not filled.
C
C     R*4   TKS1(3)    Track 1 swum parameters (phi,1/Pxy,tanlam)
C     R*4   TKS2(3)    Track 2 swum parameters
C     R*4   ERS1(6)    Track 1 swum error matrix (usual sequence)
C     R*4   ERS2(6)    Track 2 swum error matrix
C     R*4   CRXY(2)    (x,y) position of chosen crossing poing (meters)
C     R*4   ZSWUM(2)   z position of tracks 1 and 2 at crossing point (m)
C     R*4   P4PAR(4)   Parent (Px,Py,Pz,E) after swimming (GeV)
C     R*4   MPAR       Parent invariant mass (GeV)
C     R*4   PXYZ       Magnitude of parent three-momentum (GeV)
C     R*4   CHI2XY     Chi-squared (1 dof) of xy alignment of parent
C                       momentum with line from interaction point
C                       to crossing point.  See Mark III Memo
C                       referenced above.
C     R*4   CHI2Z      Chi-squared (1 dof) for the hypothesis that
C                       the two tracks have equal z coordinates at
C                       the crossing point.  **WARNING**: This
C                       quantity would have no meaning if one or
C                       both of the tracks has mzfit .ne. 3, so
C                       CHI2Z = -1.0 is returned in this case.
C     R*4   LXY        Best fit value of lab xy decay distance (m)
C     R*4   CT         c * proper decay time (m), derived from LXY,
C                       P4PAR, and the exact K0 or Lambda mass.
C
C
C                            Steve Wasserbaech  11 April 1988
C
C   In order to make this subroutine usefull at BES,we use dedx corrected
C   (Mike Kelsey's formulae) track to replace of the original dedx track
C   from trklst available at MARK III(pointed out by Bruce Lowery),and add
C   a variable  of "iversion" to deal with  the error matrix for the
C   different reconstruction version.
C
C                                         H.B.Lan
C                                         Jan. 11,1993
c   For iversion.ge.85,make use of mfit=-9 track also,pointed
c   out by Bruce Lowery
c                                         H.B.Lan
C
C   RECOVER THE GEOMETRICAL ERROR MATRIX AFTER K0 FINDING,
C                                         H.B.LAN,MAY,4,94
C
C   AS IVERSION 831 IS THE LOWER IVERSION THAN 85 AT BES,CORRECT THE
C   FORMER BUG.
C                                         H.B.LAN ,AUG.,15,1994
C-----------------------------------------------------------------------

klams0.f

C-----------------------------------------------------------------------
C
C   This file contains:
C
C    Subroutine KLAMS0(MASS1,MASS2,Q1,Q2,TK1,TK2,ER1,ER2,IRC,TKS1,
C    *TKS2,ERS1,ERS2,CRXY,ZSWUM,P4PAR,MPAR,PXYZ,CHI2XY,CHI2Z,LXY,CT)
C
C
C    KLAMS0 performs the swimming to the crossing point and other
C  related calculations.  The track parameters and error matrices
C  are input.  Only the position and size of the interaction region
C  are taken from Trklst.
C
C
C  Input:
C     R*4   MASS1,MASS2  Masses (GeV) of tracks 1 and 2
C     I*4   Q1,Q2        Charges of tracks 1 and 2
C
C                        The following track information should
C                         already be corrected for energy loss
C                         and multiple scattering:
C     R*4   TK1(6)       Track 1 parameters (phi,1/Pxy,tanlam,xm,ym,eta)
C     R*4   TK2(6)       Track 2 parameters
C     R*4   ER1(15)      Track 1 error matrix (usual sequence)
C     R*4   ER2(15)      Track 2 error matrix
C
C  Output:
C     I*4   IRC        Return code
C                       = 0  Normal completion;
C                         4  Tracks are not oppositely charged;
C                         5  Mass1, Mass2 are not consistent
C                                with pi pi or p pi;
C                         6  Parent has Pxy = 0 before swimming;
C                         7  Daughter has Pxy = 0;
C                         8  Track circles do not intersect;
C                         9  Two good crossing points found;
C                        10  Zero good crossing points found;
C                        11  Parent has Pxy = 0 after swimming;
C                        12  Arithmetic error computing uncertainty
C                                in crossing point position;
C
C                         For IRC=0 all outputs are filled.
C
C                         Fatal return codes are 1,3,5,6,7,8,11,12;
C                       none of the other outputs are filled.
C
C                         For IRC=9 the closer crossing point to the
C                       interaction point is chosen and all outputs
C                       are filled.
C
C                         For IRC=4 or 10 no swimming is done.
C                       TKS1, TKS2, ERS1, ERS2, P4PAR, MPAR, PXYZ
C                       are filled using the dE/dx corrected (unswum)
C                       track parameters.  CRXY, ZSWUM, CHI2XY,
C                       CHI2Z, LXY, CT are not filled.
C
C     R*4   TKS1(3)    Track 1 swum parameters (phi,1/Pxy,tanlam)
C     R*4   TKS2(3)    Track 2 swum parameters
C     R*4   ERS1(6)    Track 1 swum error matrix (usual sequence)
C     R*4   ERS2(6)    Track 2 swum error matrix
C     R*4   CRXY(2)    (x,y) position of chosen crossing poing (meters)
C     R*4   ZSWUM(2)   z position of tracks 1 and 2 at crossing point (m)
C     R*4   P4PAR(4)   Parent (Px,Py,Pz,E) after swimming (GeV)
C     R*4   MPAR       Parent invariant mass (GeV)
C     R*4   PXYZ       Magnitude of parent three-momentum (GeV)
C     R*4   CHI2XY     Chi-squared (1 dof) of xy alignment of parent
C                       momentum with line from interaction point
C                       to crossing point.  See Mark III Memo
C                       referenced above.
C     R*4   CHI2Z      Chi-squared (1 dof) for the hypothesis that
C                       the two tracks have equal z coordinates at
C                       the crossing point.  **WARNING**: This
C                       quantity would have no meaning if one or
C                       both of the tracks has mzfit .ne. 3, so
C                       CHI2Z = -1.0 is returned in this case.
C     R*4   LXY        Best fit value of lab xy decay distance (m)
C     R*4   CT         c * proper decay time (m), derived from LXY,
C                       P4PAR, and the exact K0 or Lambda mass.
C
C
C                            Steve Wasserbaech  11 April 1988
C
C
C  Add better protection against tangent circles.  SRW  11 Jan 1989
C
C  Add check to see of VRTSUB(3,*) is filled before using it.
C                            B. Lowery  28 Jan 1994
C  Add protection against divide by ZERO error in evaluation of 'lxy'.
C
C------------------------------------------------------------------------

klamslow.f

C-----------------------------------------------------------------------
C
C   This file contains:
C
C    Subroutine KLAMSLOW(ITK1,ITK2,MASS1,MASS2,IRC,TKS1,TKS2,ERS1,
C     *ERS2,CRXY,ZSWUM,P4PAR,MPAR,PXYZ,CHI2XY,CHI2Z,LXY,CT,IVERSION)
C
C
C    KLAMSLOW is the same as KLAMS, except that corrections for
C  dE/dx and multiple scattering are actually computed.  Use this
C  routine when the dE/dx-corrections blocks are not available in
C  Trklst track subtype 3 for the pion and/or proton hypothesis.
C
C   The heart of the calculation is placed in Subroutine Klams0,
C  allowing the user to save cpu time if the same charged tracks
C  are to be used in several combinations.  The user calls
C  SUPERSD or SUPERSDg for each track to perform the
C  dE/dx corrections, then loops over combinations calling Klams0.
C
C  Input:
C     I*4   ITK1,ITK2    Trklst track numbers
C     R*4   MASS1,MASS2  Masses (GeV) of tracks 1 and 2
C     I     IVERSION     THE RECONSTRUCTED VERSION NUMBER H.B.LAN
C  Output:
C     I*4   IRC        Return code
C                       = 0  Normal completion;
C                         1  Specified track(s) not charged;
C                         3  Specified track(s) not helix fit;
C                         4  Tracks are not oppositely charged;
C                         5  Mass1, Mass2 are not consistent
C                                with pi pi or p pi;
C                         6  Parent has Pxy = 0 before swimming;
C                         7  Daughter has Pxy = 0;
C                         8  Track circles do not intersect;
C                         9  Two good crossing points found;
C                        10  Zero good crossing points found;
C                        11  Parent has Pxy = 0 after swimming;
C                        12  Arithmetic error computing uncertainty
C                                in crossing point position;
C
C                         For IRC=0 all outputs are filled.
C
C                         Fatal return codes are 1,3,5,6,7,8,11,12;
C                       none of the other outputs are filled.
C
C                         For IRC=9 the closer crossing point to the
C                       interaction point is chosen and all outputs
C                       are filled.
C
C                         For IRC=4 or 10 no swimming is done.
C                       TKS1, TKS2, ERS1, ERS2, P4PAR, MPAR, PXYZ
C                       are filled using the dE/dx corrected (unswum)
C                       track parameters.  CRXY, ZSWUM, CHI2XY,
C                       CHI2Z, LXY, CT are not filled.
C
C     R*4   TKS1(3)    Track 1 swum parameters (phi,1/Pxy,tanlam)
C     R*4   TKS2(3)    Track 2 swum parameters
C     R*4   ERS1(6)    Track 1 swum error matrix (usual sequence)
C     R*4   ERS2(6)    Track 2 swum error matrix
C     R*4   CRXY(2)    (x,y) position of chosen crossing poing (meters)
C     R*4   ZSWUM(2)   z position of tracks 1 and 2 at crossing point (m)
C     R*4   P4PAR(4)   Parent (Px,Py,Pz,E) after swimming (GeV)
C     R*4   MPAR       Parent invariant mass (GeV)
C     R*4   PXYZ       Magnitude of parent three-momentum (GeV)
C     R*4   CHI2XY     Chi-squared (1 dof) of xy alignment of parent
C                       momentum with line from interaction point
C                       to crossing point.  See Mark III Memo
C                       referenced above.
C     R*4   CHI2Z      Chi-squared (1 dof) for the hypothesis that
C                       the two tracks have equal z coordinates at
C                       the crossing point.  **WARNING**: This
C                       quantity would have no meaning if one or
C                       both of the tracks has mzfit .ne. 3, so
C                       CHI2Z = -1.0 is returned in this case.
C     R*4   LXY        Best fit value of lab xy decay distance (m)
C     R*4   CT         c * proper decay time (m), derived from LXY,
C                       P4PAR, and the exact K0 or Lambda mass.
C
C
C                            Steve Wasserbaech  11 April 1988
C
C
C------------------------------------------------------------------------

newsquid.f

C-----------------------------------------------------------------------
C
C  This file contains:
C
C   Function   NEWSQUID(N_PART,MASS,P_SPEC,N_RES,M_RES,P_MEM,
C     *MEM,RC,CHISQ,P_FIT)
C
C
C   Newsquid is called to invoke a kinematic fit.  The total event
C energy and momentum constraints are imposed.
C
C   Before invoking Newsquid, it is necessary to call:
C      1) Psinit     to initialize
C      2) PSTRKCP, PSTRKCO, PSBMFIT,
C         Psexternal, PSDOTRKC, and/or Psdobeamfit
C                    to specify the tracks (one at a time)
C      3) PsresonanceN (N=2,3,...,10)
C                    to specify the resonances, if any (one at a time)
C
C   Do HELP MARKIII TELESIS for complete instructions.  The sequence
C described will automatically fill the input parameters of Newsquid,
C which are almost unchanged from the old Squid.  The help file also
C contains a detailed description of these old parameters.
C   For fitting examples, see TELESIS EXAMPLES.
C
C   Input:
C      I*4   N_PART       Number of particles
C      R*4   MASS(*)      Particles masses
C      I*4   P_SPEC(4,*)  See note 3 of PSTRKCP.
C      I*4   N_RES(*)     Number of resonances
C      R*4   M_RES(*)     Resonances masses
C      I*4   P_MEM(*)     Pointer into MEM for resonances
C      I*4   MEM(*)       Array for specifying tracks in resonances
C      I*4   RC           Return code from pre-squid calls
C
C  Output:
C      L*4   NEWSQUID     Return code
C                          =.TRUE.   Successful completion
C                          =.FALSE.  Fit failed--error in Newsquid
C                                     or in pre-Squid routines.
C
C      I*4   RC           Return code
C                          = 0       Successful completion
C                          = 7       Fit failed--LEQT2P found non-
C                                     positive-definite matrix or
C                                     did not converge.  (It is
C                                     normal to find this failure
C                                     once in a while.)
C                          = 8       Fit failed--chi-squared > 50.
C
C      R*4   CHISQ        chi-squared of fit
C      R*4   P_FIT(*)     Fit particle four-momenta
C
C-----------------------------------------------------------------------

psbmfit.f

C-----------------------------------------------------------------------
C
C  This file contains:
C
C   Subroutine PSBMFIT(I_SQUID,MASS_I,SPEC1,SPEC2,SPEC3,SPEC4,
C     *N_PART,MASS,P_SPEC,RC,BCHISQ,XI0)
C
C
C   PSBMFIT is called to specify one particle for a Newsquid or
C Squidjr fit.  The track parameters and error matrix (corrected
C for dE/dx and multiple scattering, and beamfit [see BMFIT MORTRAN])
C are taken from the extended Trklst track subtype 3.  The track
C information is correctly loaded into the Newsquid Common block.
C
C   The parameters of this routine are exactly the same as for
C PSTRKCP.  (See comments therein.)
C
C  Input:
C     I*4   I_SQUID      Squid track number.  See note 1 of PSTRKCP.
C     R*4   MASS_I       Particle Mass (GeV/c**2)
C     I*4   SPEC1        See note 3.
C     I*4   SPEC2        ''
C     I*4   SPEC3        ''
C     I*4   SPEC4        TRKLST track number.  See note 3.
C     I*4   N_PART       Number of particles so far
C     R*4   MASS(*)      Masses of particles
C     I*4   P_SPEC(4,*)  See note 3.
C     I*4   RC           Return code.  See note 4.
C
C Output:
C     I*4   N_PART       New number of particles
C     R*4   MASS(*)      Masses of particles
C     I*4   P_SPEC(4,*)  See note 3.
C     I*4   RC           Return code
C                         = 0       Successful completion
C                         = 1       Track was not charged!  Specify
C                                    neutral tracks with PSTRKCP
C                                    or Psexternal, and missing tracks
C                                    with PSTRKCP.
C                         = 2       Specified mass invalid; it must
C                                    correspond to e, mu, pi, K, or p,
C                                    within 5% for the info block
C                                    to be found.  mass_i is used in
C                                    the fit without adjustment.
C                         = 5       Track did not have dE/dx block
C                                    filled for the specified mass
C                                    hypothesis.
C                         = 6       Track did not have Bmfit block
C                                    filled for the specified mass
C                                    hypothesis.
C
C     R*4   BCHISQ       Bmfit chi-squared (1 degree of freedom).
C
C     R*4   XI0          SIGNED distance of closest approach (meters)
C                          to BEAM axis BEFORE the Bmfit (after dE/dx
C                          corrections).  Cut on abs(XI0) instead of RM.
C
C
C                        If (rc .ne. 0) on input, then execution
C                        is aborted and rc is left unchanged.
C
C                                    SRW  21 Apr 88
C
C-----------------------------------------------------------------------

psdobmft.f

C-----------------------------------------------------------------------
C
C  This file contains:
C
C   Subroutine PSDOBMFT(I_SQUID,MASS_I,SPEC1,SPEC2,SPEC3,SPEC4,
c     *N_PART,MASS,P_SPEC,RC,BCHISQ,XI0,IVERSION)
C
C
C   PSDOBMFT is called to specify one particle for a Newsquid
C or Squidjr fit.  The track momentum is taken from TRKLST and
C corrections for dE/dx and multiple scattering are applied
C to the track parameters and error matrix.  A beam fit is also
C performed (see BMFIT MORTRAN).  The track information is
C correctly loaded into the Squid Common block.
C
C   The parameters of this routine are exactly the same as for
C PSTRKCP.  (See comments therein.)
C
C   NOTE: In preliminary versions of Newsquid, this routine was
C called Psbeamfit.
C
C  Input:
C     I*4   I_SQUID      Squid track number.  See note 1 of PSTRKCP.
C     R*4   MASS_I       Particle Mass (GeV/c**2)
C     I*4   SPEC1        See note 3.
C     I*4   SPEC2        ''
C     I*4   SPEC3        ''
C     I*4   SPEC4        TRKLST track number.  See note 3.
C     I*4   N_PART       Number of particles so far
C     R*4   MASS(*)      Masses of particles
C     I*4   P_SPEC(4,*)  See note 3.
C     I*4   RC           Return code.  See note 4.
C
C Output:
C     I*4   N_PART       New number of particles
C     R*4   MASS(*)      Masses of particles
C     I*4   P_SPEC(4,*)  See note 3.
C     I*4   RC           Return code
C                         = 0       Successful completion
C                         = 1       Track was not charged!  Specify
C                                    neutral tracks with PSTRKCP
C                                    or Psexternal, and missing tracks
C                                    with PSTRKCP.
C                         = 3       Track did not have MFIT=2
C                         = 4       Beam fit failed.
C
C     R*4   BCHISQ       Bmfit chi-squared (1 degree of freedom).
C
C     R*4   XI0          SIGNED distance of closest approach (meters)
C                          to BEAM axis BEFORE the Bmfit (after dE/dx
C                          corrections).  Cut on abs(XI0) instead of RM.
C
C
C                        If (rc .ne. 0) on input, then execution
C                      is aborted and rc is left unchanged.
C
C                                    SRW  22 Apr 88
C
C-----------------------------------------------------------------------

psdotrkc.f

C-----------------------------------------------------------------------
C
C  This file contains:
C
C   Subroutine PSDOTRKC(I_SQUID,MASS_I,SPEC1,SPEC2,SPEC3,SPEC4,
C     *N_PART,MASS,P_SPEC,RC,IVERSION)
C
C
C   PSDOTRKC is called to specify one particle for a Newsquid
C or Squidjr fit.  The track momentum is taken from TRKLST and
C corrections for dE/dx and multiple scattering are applied
C to the track parameters and error matrix.  The track information
C is correctly loaded into the Squid Common block.
C
C   The parameters of this routine are exactly the same as for
C PSTRKCP.  (See comments therein.)
C
C   NOTE: In preliminary versions of Newsquid, this routine was
C called PSTRKCO.
C
C  Input:
C     I*4   I_SQUID      Squid track number.  See note 1 of PSTRKCP.
C     R*4   MASS_I       Particle Mass (GeV/c**2)
C     I*4   SPEC1        See note 3.
C     I*4   SPEC2        ''
C     I*4   SPEC3        ''
C     I*4   SPEC4        TRKLST track number.  See note 3.
C     I*4   N_PART       Number of particles so far
C     R*4   MASS(*)      Masses of particles
C     I*4   P_SPEC(4,*)  See note 3.
C     I*4   RC           Return code.  See note 4.
C
C Output:
C     I*4   N_PART       New number of particles
C     R*4   MASS(*)      Masses of particles
C     I*4   P_SPEC(4,*)  See note 3.
C     I*4   RC           Return code
C                         = 0       Successful completion
C                         = 1       Track was not charged!  Specify
C                                    neutral tracks with PSTRKCP
C                                    or Psexternal, and missing tracks
C                                    with PSTRKCP.
C                         = 3       Track did not have MFIT=2
C
C                        If (rc .ne. 0) on input, then execution
C                        is aborted and rc is left unchanged.
C
C                                       SRW  22 Apr 88
C
C-----------------------------------------------------------------------

psextern.f

C-----------------------------------------------------------------------
C
C  This file contains:
C
C   Subroutine PSEXTERN(I_SQUID,MASS_I,SPEC1,SPEC2,SPEC3,SPEC4,
C     *TYPE_I,TKP,ERR,  N_PART,MASS,P_SPEC,RC)
C
C
C   PSEXTERN is called to specify one particle for a Newsquid
C or Squidjr fit.  The track momentum and error matrix are provided
C as inputs to this routine.  The track information is correctly
C loaded into the New Squid common block.
C
C   When a particular track is to be used in many fits it is
C possible to save cpu seconds by calling Sundae and/or BANASPI
C once per track per mass hypothesis and saving the resulting
C track parameters.  PSEXTERN is called to load the information
C prior to each fit.
C
C   PSEXTERN is also ideal for loading the swum track parameters
C and error matrices of K0 and Lambda daughters.  The swimming
C can be performed by the Telesis routine Klams.
C
C   The parameters of this routine are DIFFERENT than those for
C PSTRKCP:
C
C  Input:
C     I*4   I_SQUID      Squid track number.  (See note 1 of PSTRKCP.)
C     R*4   MASS_I       Particle Mass (GeV/c**2)
C     I*4   SPEC1        See note 3.
C     I*4   SPEC2        ''
C     I*4   SPEC3        ''
C     I*4   SPEC4        TRKLST track number.
C     I*4   TYPE_I       Type of track:
C                          1   Charged
C                          2   Neutral (Barrel)
C                          3   Neutral (North Endcap)
C                          4   Neutral (South Endcap)
C     R*4   TKP(3)       Track parameters;  See table in note 3.
C                        For charged tracks these are:
C                               phi, 1/Pxy, tan(lam).
C     R*4   ERR(6)       Error Matrix elements.
C     I*4   N_PART       Number of particles so far
C     R*4   MASS(*)      Masses of particles
C     I*4   P_SPEC(4,*)  See note 3.
C     I*4   RC           Return code.  See note 4.
C
C Output:
C     I*4   N_PART       New number of particles
C     R*4   MASS(*)      Masses of particles
C     I*4   P_SPEC(4,*)  See note 3.
C     I*4   RC           Return code
C                         = 0       Successful completion
C
C                        If (rc .ne. 0) on input, then execution
C                        is aborted and rc is left unchanged.
C
C----------------------------------------------------------------------

psinit.f

C-----------------------------------------------------------------------
C
C  This file contains:
C
C   Subroutine PSINIT(SQRT_S,N_PART,N_RES,RC)
C
C
C   Psinit is called for initialization before a Newsquid fit.
C
C  Input:
C     R*4   SQRT_S       Center-of-Mass Energy (GeV)
C
C Output:
C     (These output parameters are really just being initialized.)
C     I*4   N_PART       Number of particles  = 0 ALWAYS;
C     I*4   N_RES        Number of resonances = 0 ALWAYS;
C     I*4   RC           Return code          = 0 ALWAYS.
C
C  Some things in the New Squid common block are initialized here too.
C
C-----------------------------------------------------------------------

psinitjr.f

C-----------------------------------------------------------------------
C
C  This file contains:
C
C   Subroutine PSINITJR(N_PART,N_RES,RC)
C
C
C   Psinitjr is called for initialization before a Squidjr fit.
C
C  Input:
C     None.
C
C Output:
C     (These output parameters are really just being initialized.)
C     I*4   N_PART       Number of particles  = 0 ALWAYS;
C     I*4   N_RES        Number of resonances = 0 ALWAYS;
C     I*4   RC           Return code          = 0 ALWAYS.
C
C  Some things in the New Squid common block are initialized here too.
C
C-----------------------------------------------------------------------

psres10.f

C-----------------------------------------------------------------------
C
C  This file contains:
C
C   Subroutine PSRES10(M_RES_I,I1,I2,I3,I4,I5,I6,I7,I8,I9,I10,
C     *  N_RES,M_RES,P_MEM,MEM,RC)
C
C
C   PSRES10 is called to specify a ten-particle resonance in
C a Newsquid or Squidjr fit.  I1, I2,...,I10 are the Squid particle
C numbers of the particles in the resonance.  See PSRES2 for
C an explanation of the other arguments.
C
C-----------------------------------------------------------------------

psres2.f

C-----------------------------------------------------------------------
C
C  This file contains:
C
C   Subroutine PSRES2
C
C
C   PSRES2 is called to specify a two-particle resonance in
C a Newsquid or Squidjr fit.  N_RES, M_RES, P_MEM, and MEM are
C filled by this routine.  N_RES should have been initialized by
C Psinit or Psinitjr.  The others need no initialization.
C Do HELP MARKIII TELESIS to see the old Squid documentation which
C explains N_RES, M_RES, P_MEM, etc.
C
C  Input:
C     I*4   M_RES_I    Resonance mass.
C     I*4   I1         Squid particle number of first particle
C     I*4   I2         Squid particle number of second particle.
C     I*4   N_RES      Number of resonances
C     I*4   M_RES(*)   Masses of resonances
C     I*4   P_MEM(*)   Pointer into MEM for resonances
C     I*4   MEM(*)     Array for specifying tracks in resonances
C     I*4   RC         Return code
C
C
C Output:
C     I*4   N_RES      New number of resonances
C     I*4   M_RES(*)   Masses of resonances
C     I*4   P_MEM(*)   Pointer into MEM for resonances
C     I*4   MEM(*)     Array for specifying tracks in resonances
C     I*4   RC         Return code
C
C  If rc=0 is input, PSRES2 will always give rc=0 for output.
C  (This routine has no non-fatal failures.)
C
C  Use PsresonanceN for N-particle resonances (N=2,3,...,10).
C  The other Psresonance routines are identical to this one, except
C  for the number of particles being specified.
C
C-----------------------------------------------------------------------

psres3.f

C-----------------------------------------------------------------------
C
C  This file contains:
C
C   Subroutine PSRES3
C
C
C   PSRES3 is called to specify a three-particle resonance in
C a Newsquid or Squidjr fit.  I1, I2, and I3 are the Squid particle
C numbers of the particles in the resonance.  See PSRES2 for
C an explanation of the other arguments.
C
C-----------------------------------------------------------------------

psres4.f

C-----------------------------------------------------------------------
C
C  This file contains:
C
C   Subroutine PSRES4
C
C
C   PSRES4 is called to specify a four-particle resonance in
C a Newsquid or Squidjr fit.  I1, I2, I3, and I4 are the Squid particle
C numbers of the particles in the resonance.  See PSRES2 for
C an explanation of the other arguments.
C
C-----------------------------------------------------------------------

psres5.f

C-----------------------------------------------------------------------
C
C  This file contains:
C
C   Subroutine PSRES5
C
C
C   PSRES5 is called to specify a five-particle resonance in
C a Newsquid or Squidjr fit.  I1, I2,...,I5 are the Squid particle
C numbers of the particles in the resonance.  See PSRES2 for
C an explanation of the other arguments.
C
C-----------------------------------------------------------------------

psres6.f

C-----------------------------------------------------------------------
C
C  This file contains:
C
C   Subroutine PSRES6
C
C
C   PSRES6 is called to specify a six-particle resonance in
C a Newsquid or Squidjr fit.  I1, I2,...,I6 are the Squid particle
C numbers of the particles in the resonance.  See PSRES2 for
C an explanation of the other arguments.
C
C-----------------------------------------------------------------------

psres7.f

C-----------------------------------------------------------------------
C
C  This file contains:
C
C   Subroutine PSRES7
C
C
C   PSRES7 is called to specify a seven-particle resonance in
C a Newsquid or Squidjr fit.  I1, I2,...,I7 are the Squid particle
C numbers of the particles in the resonance.  See PSRES2 for
C an explanation of the other arguments.
C
C-----------------------------------------------------------------------

psres8.f

C-----------------------------------------------------------------------
C
C  This file contains:
C
C   Subroutine PSRES8
C
C
C   PSRES8 is called to specify an eight-particle resonance in
C a Newsquid or Squidjr fit.  I1, I2,...,I8 are the Squid particle
C numbers of the particles in the resonance.  See PSRES2 for
C an explanation of the other arguments.
C
C-----------------------------------------------------------------------

psres9.f

C-----------------------------------------------------------------------
C
C  This file contains:
C
C   Subroutine PSRES9
C
C
C   PSRES9 is called to specify a nine-particle resonance in
C a Newsquid or Squidjr fit.  I1, I2,...,I9 are the Squid particle
C numbers of the particles in the resonance.  See PSRES2 for
C an explanation of the other arguments.
C
C-----------------------------------------------------------------------

pstrkco.f

C-----------------------------------------------------------------------
C
C  This file contains:
C
C   Subroutine PSTRKCO
C
C
C   PSTRKCO is called to specify one particle for a Newsquid or
C Squidjr fit.  The track parameters and error matrix (corrected
C for dE/dx and multiple scattering) are taken from the extended
C Trklst track subtype 3.  The track information is correctly
C loaded into the Newsquid Common block.
C
C   The parameters of this routine are exactly the same as for
C PSTRKCP.  (See comments therein.)
C
C  Input:
C     I*4   I_SQUID      Squid track number.  See note 1 of PSTRKCP.
C     R*4   MASS_I       Particle Mass (GeV/c**2)
C     I*4   SPEC1        See note 3.
C     I*4   SPEC2        ''
C     I*4   SPEC3        ''
C     I*4   SPEC4        TRKLST track number.  See note 3.
C     I*4   N_PART       Number of particles so far
C     R*4   MASS(*)      Masses of particles
C     I*4   P_SPEC(4,*)  See note 3.
C     I*4   RC           Return code.  See note 4.
C
C Output:
C     I*4   N_PART       New number of particles
C     R*4   MASS(*)      Masses of particles
C     I*4   P_SPEC(4,*)  See note 3.
C     I*4   RC           Return code
C                         = 0       Successful completion.
C                         = 1       Track was not charged!  Specify
C                                    neutral tracks with PSTRKCP
C                                    or Psexternal, and missing tracks
C                                    with PSTRKCP.
C                         = 2       Specified mass invalid; it must
C                                    correspond to e, mu, pi, K, or p,
C                                    within 5% for the info block
C                                    to be found.  mass_i is used in
C                                    the fit without adjustment.
C                         = 5       Track did not have dE/dx block
C                                    filled for the specified mass
C                                    hypothesis.
C
C                        If (rc .ne. 0) on input, then execution
C                        is aborted and rc is left unchanged.
C
C-----------------------------------------------------------------------

pstrkcp.f

C-----------------------------------------------------------------------
C
C  This file contains:
C
C   Subroutine PSTRKCP
C
C
C   PSTRKCP is called to specify one particle for a Newsquid or
C Squidjr fit.  The track momentum is taken directly from TRKLST.
C No dE/dx corrections are applied.  The track information is
C correctly loaded into the New Squid common block.
C
C   Use this routine for specifying neutral tracks and missing tracks.
C
C   When this routine is called for an appropriate track, N_PART
C is incremented and MASS() and P_SPEC() are filled with the
C track information.  The track parameters are placed in New Squid
C common. See Squid Mortran G to see the old Squid documentation which
C explains N_PART, MASS, P_SPEC, etc.
C
C  Input:
C     I*4   I_SQUID      Squid particle number.  See note 1 below.
C     R*4   MASS_I       Particle Mass (GeV/c**2).  See note 2.
C     I*4   SPEC1        See note 3.
C     I*4   SPEC2        ''
C     I*4   SPEC3        ''
C     I*4   SPEC4        TRKLST track number.  See note 3.
C     I*4   N_PART       Number of particles so far
C     R*4   MASS(*)      Masses of particles
C     I*4   P_SPEC(4,*)  See note 3.
C     I*4   RC           Return code.  See note 4.
C
C Output:
C     I*4   N_PART       New number of particles
C     R*4   MASS(*)      Masses of particles
C     I*4   P_SPEC(4,*)  See note 3.
C     I*4   RC           Return code
C
C Notes:
C
C 1.  When a set of tracks is being specified for a Squid fit, the
C   first call to PSTRKCP, PSTRKCO, PSBMFIT, PSEXTERN,
C   psdoTRKcorr, or PSDOBMFT must have I_SQUID=1, the second
C   must have I_SQUID=2, etc. This safety feature prevents
C   certain bugs from slipping through undetected.  It also
C   makes your resonance definitions easy to trace.
C
C 2.  It is possible to impose an equal mass constraint in Newsquid
C   (but not in Squidjr).  There are two cases:
C          a)  M(resonance1) =  M(resonance2)
C       or b)  M(resonance)  =  M(missing particle)
C   The mass to which the two object converge is an output of the
C   fit.  To impose such a constraint, specify the mass of the
C   objects to be -1.0 GeV.  (The missing particle mass is specified
C   in PSTRKCP, and the resonance mass in one of the Psresonance
C   routines.)
C
C 3.  In the old SQUID, the array PSPEC allowed the user to specify
C   missing track parameters (e.g. the energy of a K-long which
C   produced a shower in the calorimeter).  Each track had four
C   PSPEC quantities.  The first three corresponded to the three
C   track parameters used in Squid.  The fourth was the TRKLST
C   track number, or zero for a missing track.
C
C                                        ---Track parameters ---
C                              type      1          2          3
C
C     Charged track              1      phi       1/Pxy     tan(lam)
C     Neutral (Barrel)           2      phi      tan(lam)   sqrt(E)
C     Neutral (North Endcap)     3      x/z        y/z      sqrt(E)
C     Neutral (South Endcap)     4      x/z        y/z      sqrt(E)
C
C   [Type is an argument of PSEXTERN.]
C   Actually the only parameter allowed to be missing in Newsquid
C   (as in the old Squid) is sqrt(E) for neutrals.
C
C   When specifying a track in Newsquid one places the PSPEC
C   quantities in Spec1 through Spec4.  These values are
C   automatically copied into the appropriate slots of P_SPEC
C   in preparation for the fit.
C
C 4.  The return codes of the pre-fitting and fitting routines
C   are designed to allow uncluttered passages of fitting code.
C   [Pre-fitting routines are those having names beginning with
C   `PS' (= Pre-Squid).]
C
C   It is NOT necessary to check for a successful return from
C   each call.  The rc variable is zeroed by psinit or psinitjr.
C   Any of the subsequent pre-fitting or fitting routines will
C   return immediately if they receive a non-zero rc as input.
C   (The value of rc will not be changed.)  If any failure occurs
C   in the pre-squid calls, the fitting function subprogram
C   (Newsquid or Squidjr) will be returned as .false..
C   See TELESIS EXAMPLES for examples.
C
C   If rc=0 is input, PSTRKCP will always give rc=0 for output.
C   (This routine has no non-fatal failures.)  Fatal errors are
C   described in Subroutine SQUIDFAT.
C
C-----------------------------------------------------------------------
C
C     Add x0, y0 and z0 etc. coords. of center of interaction region
C   to neutral tracks. The coords. can be extracted from Subtype 8 or
C   from the common block RUNPRM which was filled by calling subroutine
C   BMSTAT. Here I use commom block RUNPRM.
C
C                         Jan. 7, 1993,  Shenjian Chen
C     For iversion less than 8,use Mike Kelsey's dedx correction to
c    modify the track parameters
c
c                         Jan.11,1994
c                                         Huibin Lan
C-----------------------------------------------------------------------

squidfat.f

C-----------------------------------------------------------------------
C
C  This file contains:
C
C   Subroutine SQUIDFAT
C
C
C  SQUIDFAT (an internal Squid routine) is called when a
C  fatal error is detected.  A message is printed and Eojob is
C  called to terminate execution.
C
C  Input:
C     I*4   CODE       Error number
C
C     (The definitions of the following inputs depends on the
C     error code being handled.)
C     I*4   I1
C     I*4   I2
C
C
C  SRW  6 Jul 88  Changed evsums call to dstprt.
C
C----------------------------------------------------------------------

squidjr.f

C-----------------------------------------------------------------------
C
C  This file contains:
C
C   Function   squidjr
C
C
C   SQUIDJR is called to invoke a kinematic fit.  The total event
C energy and momentum constraints are NOT imposed!  Only resonance
C constraints are allowed.
C
C   Before invoking Squidjr, it is necessary to call:
C      1) Psinitjr   to initialize
C      2) PSTRKCP, PSTRKCO, PSBMFIT,
C         Psexternal, PSDOTRKC, and/or Psdobeamfit
C                    to specify the tracks (one at a time)
C      3) PsresonanceN (N=2,3,...,10)
C                    to specify the resonances (one at a time)
C
C   Do HELP MARKIII TELESIS for complete instructions.  The sequence
C described will automatically fill the input parameters of Squidjr,
C which are almost unchanged from the old Squid.  The help file also
C contains a detailed description of these old parameters.
C   For fitting examples, see TELESIS EXAMPLES.
C
C   Input:
C      I*4   N_PART       Number of particles
C      R*4   MASS(*)      Particles masses
C      I*4   P_SPEC(4,*)  See note 3 of PSTRKCP.
C      I*4   N_RES(*)     Number of resonances
C      R*4   M_RES(*)     Resonances masses
C      I*4   P_MEM(*)
C      I*4   MEM(*)
C      I*4   RC           Return code
C
C  Output:
C      I*4   RC           Return code
C                          = 0       Successful completion
C                          = 7       Fit failed--LEQT2P found non-
C                                      positive-definite matrix or
C                                      did not converge.  (It is
C                                      normal to find this failure
C                                      once in a while.)
C                          = 8       Fit failed--chi-squared > 50.
C
C      R*4   CHISQ        chi-squared of fit
C      R*4   P_FIT(*)     Fit particle four-momenta
C
C
C Things to remember:
C
C 1.  The cms energy is not used in Squidjr.
C
C 2.  Missing tracks and showers with missing sqrt(E) are not allowed.
C
C 3.  Equal mass constraints are not allowed.
C
C 4.  The number of degrees of freedom of CHISQ is just N_RES
C   (= number of constraints).
C
C-----------------------------------------------------------------------

squidpul.f

C-----------------------------------------------------------------------
C
C  This file contains:
C
C   Subroutine SQUIDPUL
C
C
C   SQUIDPUL is used to calculate the pulls from a successful
C Newsquid or Squidjr fit.  Relevant formulae appear in Frodesen
C et al., `Probability and Statistics in Particle Physics,' pp.
C 283, 289, and 290.
C
C
C  Input:
C     I*4   P_SPEC(4,*)  See note 3 of PSTRKCP.
C     I*4   RC           Return code for pre-squid calls and Newsquid
C                         or Squidjr.  If RC .ne. 0 on input, the
C                         execution of this routine will be aborted
C                         and PULL(i) will be set to 999 for i=1,NVAR.
C                         (See note 1 below.)
C
C Output:
C     I*4   RC           Return code
C                         = 0       Successful completion
C                         = 9       LEQT2P failed to compute
C                                    new error matrix.  NOTICE:
C                                    This does not mean that the
C                                    fit was bad.  Check the value
C                                    of Newsquid or Squidjr to
C                                   determine whether the fit was
C                                   successful.
C
C     R*4   PULL(*)      Array of pulls.  The dimension of this
C                        array should be at least NVAR; there will be
C                        pulls placed in cells 1,2,3,...,NVAR.
C                        See note 1.
C
C Notes:
C 1.  Normally NVAR=3*N, where N is the number of MEASURED tracks
C     [If  you have omitted a neutral shower energy measurement
C     via SPEC3=0, there will be 3*N - 1 pulls.]
C
C     The sequence of pulls in the array will be:
C
C            |  Track 1  |  Track 2  | . . . |    Track N     |
C     Index:   1   2   3   4   5   6   . . .   3N-2  3N-1  3N
C
C     [If you have used SPEC3=3, the track with the missing energy
C     will have only 2 pulls.  If you have specified a completely
C     missing track, that track will be skipped completely in the
C     array of pulls.]
C
C                                  SRW  22 Apr 88
C
C-----------------------------------------------------------------------

sundae.f

C----------------------------------------------------------------------
C
C  This file contains:
C
C     Subroutine Sundae
C
C
C    Sundae is used to recover the uncorrected
C  error matrix for a charged track (ICECREAM) and apply
C  corrections to the track parameters and error matrix
C  corresponding to a particular mass hypothesis (FUDGE).
C
C     Input:
C        I*4  ITRACK      TRKLST track number
C        R*4  XMAS        desired mass hypothesis
C        I    IVERSION    THE RECONSTRUCTED VERSION ! LANHB
C     Output:
C        I*4  IRC         Return code:
C                          0  Normal completion
C                          1  Track is not charged!
C                          2  (Applies to SUNDAEG.)
C                          3  Track does not have MFIT=2.
C
C        I*4  ICHG        Track charge
C        R*4  TKP(3)      phi, 1/Pxy, tan(lambda)
C        R*4  ERR(6)      1st 6 elements of error matrix
C
C
C       If the geometric (undamaged) error matrix is available
C     in your data set, you may want to use the subroutine
C     SUNDAEG.  This is faster and gives more accurate results
C     since the error matrix does not need to undergo the
C     approximate ICECREAM process.
C
C       See SUBROUTINE SUNDAE0 below if you want corrected track
C     momenta but don't need the error matrices.
C
C                                   SRW  20 Apr 88
C
C       FOR IVERSION.GE.85,WE MAKE USE OF MFIT=-19 AND MFIT=-9 TRACK
C     ALSO.                         H.B.LAN  FEB.23,1994
C
C   AS IVERSION 831 IS THE LOWER IVERSION THAN 85 AT BES,CORRECT THE
C   FORMER BUG.
C                                         H.B.LAN ,AUG.,15,1994
C----------------------------------------------------------------------

sundae0.f

C----------------------------------------------------------------------
C
C  This file contains:
C
C    Subroutine Sundae0
C
C
C   Sundae0 is used to obtain charged track parameters
C   corrected for dE/dx; no error matrix corrections are done.
C
C     Input:
C        I*4  ITRACK      TRKLST track number
C        R*4  XMAS        desired mass hypothesis
C
C     Output:
C        I*4  IRC         Return code:
C                          0  Normal completion
C                          1  Track is not charged!
C        I*4  ICHG        Track charge
C        R*4  TKP(3)      phi, 1/Pxy, tan(lam)
C
C                                    20 Apr 88  SRW
C
C----------------------------------------------------------------------

sundaeg.f

C----------------------------------------------------------------------
C
C  This file contains:
C
C     Subroutine Sundaeg
C
C
C    Sundaeg is used to read the geometric error
C  matrix for a charged track from TRKLST and apply
C  corrections to the track parameters and error matrix
C  corresponding to a particular mass hypothesis (FUDGE).
C
C     Input:
C        I*4  ITRACK      TRKLST track number
C        R*4  XMAS        desired mass hypothesis
C        I    IVERSION    THE RECONSTRUCTED VERSION  ! LANHB
C     Output:
C        I*4  IRC         Return code:
C                          0  Normal completion
C                          1  Track is not charged!
C                          2  Geometric error matrix not found.
C                          3  Track does not have MFIT=2.
C
C        I*4  ICHG        Track charge
C        R*4  TKP(3)      phi, 1/Pxy, tan(lam)
C        R*4  ERR(6)      1st 6 elements of error matrix
C
C                                    SRW 20 Apr 88
C       FOR IVERSION.GE.85,WE MAKE USE OF MFIT=-19 AND MFIT=-9 TRACK
C     ALSO.                         H.B.LAN  FEB.23,1994
C   AS IVERSION 831 IS THE LOWER IVERSION THAN 85 AT BES,CORRECT THE
C   FORMER BUG.
C                                         H.B.LAN ,AUG.,15,1994
C----------------------------------------------------------------------

superfug.f

C-----------------------------------------------------------------------
C
C  This file contains:
C
C   Subroutine SUPERFUG
C
C
C   SUPERFUG is used to apply dE/dx corrections to a charged
C  track's momentum, and to correct the track error matrix for
C  dE/dx and multiple scattering, according to a specified
C  mass hypothesis.
C
C  Input:
C     R*4   XMAS        Mass hypothesis (GeV/c**2);
C     R*4   TKPAR(6)    Vector of track parameters:
C                         phi, k (=1/Pxy), s (=tanlam), xm, ym, zm;
C     I*4   ICHG        Charge;
C     R*4   BC          = 0.029979 * B field at (xm,ym,zm);
C     R*4   PM(5)       Vector containing Px, Py, Pz, Pxyz, Pxy;
C     R*4   SECLM(4)    Ith element = sec(lambda)**I;
C     R*4   ERROR(15)   Error matrix (taken from TRK(I3+20) et seq.);
C     R*4   RM          = sqrt(xm**2 + ym**2);
C     I*4   MFIT        = Least-significant byte of MFITWD;
C     I*4   MZFIT       = Next-to-least-significant byte of MFITWD;
C     I*4   IFIT(2)     R-phi layer maps = ITRK(I3+35), ITRK(I3+36);
C     I*4   CD          CD CODE = ITRK(I3+40).
C
C Output:
C           TKPAR(6)    Corrected track parameters;
C           ERROR(15)   Corrected error matrix.
C
C   The following track parameters are corrected: phi, k, xm, and ym.
C
C--------------------------------------------------------------------C
C
C  Miniature set of COULOM-like constants:
C
C  U        = exponent in D. Hitlin formula
C  XMCS     = (0.0141 GeV)**2 for Multiple Scattering
C  NSCATR   = number of scatterers under consideration
C
C               In the following arrays, the scatterers
C               are numbered as follows:
C                  1. Beampipe
C                  2. Little Chamber
C                  3. Inner wall of Main DC
C                  4. Gas and wires of DC layers 2-5
C
C  RSCATR   = radius of mid-scatterer from beam axis
C  DPMIN    = d(momentum) for each scatterer (from Hitlin)
C               for a minimum-ionizing particle
C  ALPHAT   = total alpha*t for each scatterer, determined
C               from COULOM COMMON.  These numbers are
C               not necessarily consistent with the DPMIN
C               values, but it's not worth worrying about.
C  RADL     = number of radiation lengths in each scatterer
C  A seperate subroutine at BES          H.B.Lan,  Jan. 14,1994
C
C  B. Lowery Cosmetic fix to range error when INDEX1 = 940
C            just expand dimension of f1,f2 by 1.
C            May 18, 1994
C
C--------------------------------------------------------------------C

supergas.f

C--------------------------------------------------------------------C
C
C  This file contains:
C
C   Subroutine SUPERGAS
C   Block Data Gscatrcoulom
C
C
C     Routine to perform corrections to the charged track
C     error matrix ERROR to account for multiple scattering
C     in the drift chamber gas. The technique is based on
C     R.L. Gluckstern's article: NIM 24 (1963) pp.381-389.
C
C     Jan. 5, 1985  G. Blaylock
C
C
C  Input:
C     R*4   XMAS        Mass hypothesis (GeV/c**2);
C     R*4   TKPAR(6)    Vector of track parameters:
C                         phi, k (=1/Pxy), s (=tanlam), xm, ym, zm;
C     I*4   ICHG        Charge;
C     R*4   BC          = 0.029979 * B field at (xm,ym,zm);
C     R*4   PM(5)       Vector containing Px, Py, Pz, Pxyz, Pxy;
C     R*4   SECLM(4)    Ith element = sec(lambda)**I;
C     R*4   ERROR(15)   Error matrix (taken from TRK(I3+20) et seq.);
C     R*4   RM          = sqrt(xm**2 + ym**2);
C     I*4   MFIT        = Least-significant byte of MFITWD;
C     I*4   MZFIT       = Next-to-least-significant byte of MFITWD;
C     I*4   IFIT(2)     R-phi layer maps = ITRK(I3+35), ITRK(I3+36);
C     I*4   CD          CD CODE = ITRK(I3+40).
C
C  Output:
C           ERROR(15)   Error matrix, corrected for gas scattering.
C
C
C     Extend to full error matrix.
C             1987  SRW
C
C     Bug fix.  Statement label 410 was in the wrong place.
C             17 Feb 88  SRW
C
C     Fix hexadecimal constant initializer
C                                 B. Lowery  Nov 10, 1993
C
C     93/12/16 J. Panetta
C       Changed Include statements to ALL CAPS to conform with standard
C       (and not to confuse converters)
C
C     93/12/??  Unspecified changes by Lan Huibin   (comment added by Panetta)
C
C     94/01/04 J. Panetta
C       Performed EXACTLY THE SAME CHANGES I did on 93/12/16 because
C       Lan Huibin didn't include them when he changed it earlier.
C
C--------------------------------------------------------------------C

superice.f

C---------------------------------------------------------------
C
C  This file contains:
C
C   Subroutine SUPERICE
C   Subroutine positdef
C
C
C      SUPERICE is used to recover the geometric error
C    matrix, starting from the information given in TRKLST
C    Track Subtype 3 (Version 3).  Starting with Version 4
C    of the reconstruction code, the geometric error matrix
C    is SAVED, so this routine will be unnecessary.
C
C      The routine Waiter is designed to look up all of the
C    information required for input to this routine.
C
C      This is an extended version of ICECREAM, again based
C    on G. Blaylock's UNCRAP.  An attempt is made to repair
C    ALL of the damage done by ISCATR and PQSCTT.
C
C  Input:
C     R*4   TKPAR(6)    Vector of track parameters:
C                         phi, k (=1/Pxy), s (=tanlam), xm, ym, zm;
C     I*4   ICHG        Charge;
C     R*4   BC          = 0.029979 * B field at (xm,ym,zm);
C     R*4   PM(5)       Vector containing Px, Py, Pz, Pxyz, Pxy;
C     R*4   SECLM(4)    Ith element = sec(lambda)**I;
C     R*4   ERROR(15)   Error matrix (taken from TRK(I3+20) et seq.);
C     R*4   RM          = sqrt(xm**2 + ym**2);
C     I*4   MFIT        = Least-significant byte of MFITWD;
C     I*4   MZFIT       = Next-to-least-significant byte of MFITWD;
C     I*4   IFIT(2)     R-phi layer maps = ITRK(I3+35), ITRK(I3+36)!MIII;
C                       THIS IS NOT TRUE FOR BES,FOR BES,
C                       ITRK(I3+35), ITRK(I3+36)
C     I*4   CD          CD CODE = ITRK(I3+40).
C     I     IVERSION    THE RUNCONSTRUCTED VERSION   ! LANHB
C  Output:
C           ERROR(15)   Estimate of geometric error matrix.
C
C
C                                  SRW  14 Sep 1987
C
C      Guarantee error matrix positive-definiteness.
C
C                                  SRW   6 Oct 1987
C---------------------------------------------------------------
C
C  Try to adjust the subroutine to the BES enviroment
C
C                                  H.B.LAN  20 AUG 1992
C
C  AS THE GEOMETRIC ERROR MATRIX IS AVAILABLE FROM
C  IVERSION.GE.8,SO THIS SUBROUTINE IS UNNECESSARY FOR
C  IVERSION.GE.8!
C                                  H.B.LAN  29,JAN,1994
C
C
C---------------------------------------------------------------

supersd.f

C----------------------------------------------------------------------
C
C  This file contains:
C
C    Subroutine SUPERSD
C
C
C  SUPERSD is just like Sundae, except the full
C  TKPAR and ERR vectors are returned.
C
C     Input:
C        I*4  ITRACK      TRKLST track number
C        R*4  XMAS        desired mass hypothesis
C        I    IVERSION    THE RECONSTRUCTED VERSION   ! LANHB
C     Output:
C        I*4  IRC         Return code:
C                          0  Normal completion
C                          1  Track is not charged!
C                          2  (Applies to SUPERSDG.)
C                          3  Track does not have MFIT=2.
C
C        I*4  ICHG        Charge
C        R*4  TKPAR(6)    phi,k,s,xm,ym,eta
C        R*4  ERROR(15)   error matrix
C
C                                    SRW  7 Oct 87
C
C       FOR IVERSION.GE.85,WE MAKE USE OF MFIT=-19 AND MFIT=-9 TRACK
C     ALSO.                         H.B.LAN  FEB.23,1994
C   AS IVERSION 831 IS THE LOWER IVERSION THAN 85 AT BES,CORRECT THE
C   FORMER BUG.
C                                         H.B.LAN ,AUG.,15,1994
C----------------------------------------------------------------------

supersdg.f

C----------------------------------------------------------------------
C
C  This file contains:
C
C    Subroutine SUPERSDG
C
C
C  SUPERSDG is used to read the geometric error
C  matrix for a charged track from TRKLST and apply
C  corrections to the track parameters and error matrix
C  corresponding to a particular mass hypothesis (FUDGE).
C
C     Input:
C        I*4  ITRACK      TRKLST track number
C        R*4  XMAS        desired mass hypothesis
C        I                THE RECONSTRUCTED VERSION    ! LANHB
C     Output:
C        I*4  IRC         Return code:
C                          0  Normal completion
C                          1  Track is not charged!
C                          2  Geometric error matrix not found.
C                          3  Track does not have MFIT=2.
C
C        I*4  ICHG        Charge
C        R*4  TKPAR(6)    phi,k,s,xm,ym,eta
C        R*4  ERROR(15)   error matrix
C
C                                    SRW  7 Oct 87
C
C
C       FOR IVERSION.GE.85,WE MAKE USE OF MFIT=-19 AND MFIT=-9 TRACK
C     ALSO.                         H.B.LAN  FEB.23,1994
C   AS IVERSION 831 IS THE LOWER IVERSION THAN 85 AT BES,CORRECT THE
C   FORMER BUG.
C                                         H.B.LAN ,AUG.,15,1994
C----------------------------------------------------------------------

telever.f

C----------------------------------------------------------------------
C
C  This file contains:
C
C   Subroutine TELEVER
C   Block Data Firstcall
C
C
C       TELEVER print a message the first time it is called
C     in a job.
C
C                                   SRW  23 Apr 88
C
C----------------------------------------------------------------------

waiter.f

C------------------------------------------------------------------------------
C
C  This file contains:
C
C   Subroutine Waiter
C
C    WAITER must be called before ordering a SUNDAE.  This routine
C  looks up all of the TRKLST info required for correcting a track's
C  parameters and ERROR matrix for the effects of energy loss and
C  multiple scattering.  Some frequently-used quantities are also
C  calculated.
C
C  Input:
C     I*4   ITRACK      TRKLST track number
C
C  Output:
C     I*4   IRC         Return code:
C                        0  Normal completion
C                        1  ITRACK is not a charged track.
C                        2  See discussion of WAITERG below.
C
C     R*4   TKPAR(6)    Vector of track parameters:
C                         phi, k (=1/Pxy), s (=tanlam), xm, ym, zm;
C     I*4   ICHG        Charge;
C     R*4   BC          = 0.029979 * B field at (xm,ym,zm);
C     R*4   PM(5)       Vector containing Px, Py, Pz, Pxyz, Pxy;
C     R*4   SECLM(4)    Ith element = sec(lambda)**I;
C     R*4   ERROR(15)  ERROR matrix (taken from TRK(I3+20) et seq.);
C     R*4   RM          = sqrt(xm**2 + ym**2);
C     I*4   MFIT        = Least-significant byte of MFITWD;
C     I*4   MZFIT       = Next-to-least-significant byte of MFITWD;
C     I*4   IFIT(2)     R-phi layer maps = ITRK(I3+35), ITRK(I3+36);
CLHB        IFIT(2)     ITRK(I3+37), ITRK(I3+38)FOR BES
C     I*4   CD          CD CODE = ITRK(I3+40).
C
C
C    The Subroutine WAITERG may be used when the geometric ERROR
C  matrix is available in Track Subtype 3 (starting with Version 4
C  of the reconstruction code).  Usually this is preferable because
C  it allows faster and more accurate corrections to be made.  The
C  parameters for WAITERG are identical to those of WAITER except
C  ERROR (the ERROR matrix) is taken from TRK(I3+44) et seq.
C  If the routine can't find the geometric ERROR matrix, you will
C  get IRC = 2.
C
C    The Subroutine WAITER0 looks up the information needed for
C  track parameter corrections only.  Again the track number ITRACK
C  is input, and IRC, TKPAR, ICHG, BC, PM, and SECLM (as defined above)
C  are output.
C
C                                      23 Sep 87  SRW
C
C  Move WAITERG, WAITER0 into separate routines.
C                                       Lan Huibin, B. Lowery
C                                       Jan 14, 1994
C
C  FIX ONE BUG RELATED TO MFIT=-19,AND 9,POINTED OUT BY LIXN AND B.LOWERY
C                                       HUIBIN LAN MAY,26,1994
C------------------------------------------------------------------------------

waiter0.f

C------------------------------------------------------------------------------
C
C  This file contains:
C
C   Subroutine Waiter0
C
C    This subroutine is related to routine WAITER which looks up all of
C  the TRKLST info required for correcting a track's parameters and ERROR
C  matrix for the effects of energy loss and multiple scattering.  Waiter
C  also computes some frequently used quantities.
C
C    Subroutine WAITER0 looks up the information needed for
C  track parameter corrections only.  Again the track number ITRACK
C  is input, and IRC, TKPAR, ICHG, BC, PM, and SECLM (as defined above)
C  are output.
C
C    Subroutine WAITERG may be used when the geometric ERROR
C  matrix is available in Track Subtype 3 (starting with Version 4
C  of the reconstruction code).  Usually this is preferable because
C  it allows faster and more accurate corrections to be made.  The
C  parameters for WAITERG are identical to those of WAITER except
C  ERROR (the ERROR matrix) is taken from TRK(I3+44) et seq.
C  If the routine can't find the geometric ERROR matrix, you will
C  get IRC = 2.
C
C                                      23 Sep 87  SRW
C
C  Input:
C     I*4   ITRACK      TRKLST track number
C
C  Output:
C     I*4   IRC         Return code:
C                        0  Normal completion
C                        1  ITRACK is not a charged track.
C                        2  See discussion of WAITERG below.
C
C     R*4   TKPAR(6)    Vector of track parameters:
C                         phi, k (=1/Pxy), s (=tanlam), xm, ym, zm;
C     I*4   ICHG        Charge;
C     R*4   BC          = 0.029979 * B field at (xm,ym,zm);
C     R*4   PM(5)       Vector containing Px, Py, Pz, Pxyz, Pxy;
C     R*4   SECLM(4)    Ith element = sec(lambda)**I;
C     R*4   ERROR(15)  ERROR matrix (taken from TRK(I3+20) et seq.);
C     R*4   RM          = sqrt(xm**2 + ym**2);
C     I*4   MFIT        = Least-significant byte of MFITWD;
C     I*4   MZFIT       = Next-to-least-significant byte of MFITWD;
C     I*4   IFIT(2)     R-phi layer maps = ITRK(I3+35), ITRK(I3+36);
CLHB        IFIT(2)     ITRK(I3+37), ITRK(I3+38)FOR BES
C     I*4   CD          CD CODE = ITRK(I3+40).
C
C
C  Modifications:
C  ---------------------------------------------------------------------------
C  ??????????? Lan Huibin    Modified for use in BES
C
C  Jan 13,1994 Bruce Lowery  Moved WAITER, WAITERG, WAITER0, and DCRHOL
C                            into separate files.
C
C
C------------------------------------------------------------------------------

waiterg.f

C------------------------------------------------------------------------------
C
C  This file contains:
C
C   Subroutine Waiterg
C
C    This subroutine is related to routine WAITER which looks up all of
C  the TRKLST info required for correcting a track's parameters and ERROR
C  matrix for the effects of energy loss and multiple scattering.  Waiter
C  also computes some frequently used quantities.
C
C    Subroutine WAITERG may be used when the geometric ERROR
C  matrix is available in Track Subtype 3 (starting with Version 4
C  of the reconstruction code).  Usually this is preferable because
C  it allows faster and more accurate corrections to be made.  The
C  parameters for WAITERG are identical to those of WAITER except
C  ERROR (the ERROR matrix) is taken from TRK(I3+44) et seq.
C  If the routine can't find the geometric ERROR matrix, you will
C  get IRC = 2.
C
C    Subroutine WAITER0 looks up the information needed for
C  track parameter corrections only.  Again the track number ITRACK
C  is input, and IRC, TKPAR, ICHG, BC, PM, and SECLM (as defined above)
C  are output.
C                                      23 Sep 87  SRW
C
C  Input:
C     I*4   ITRACK      TRKLST track number
C
C  Output:
C     I*4   IRC         Return code:
C                        0  Normal completion
C                        1  ITRACK is not a charged track.
C                        2  See discussion of WAITERG below.
C
C     R*4   TKPAR(6)    Vector of track parameters:
C                         phi, k (=1/Pxy), s (=tanlam), xm, ym, zm;
C     I*4   ICHG        Charge;
C     R*4   BC          = 0.029979 * B field at (xm,ym,zm);
C     R*4   PM(5)       Vector containing Px, Py, Pz, Pxyz, Pxy;
C     R*4   SECLM(4)    Ith element = sec(lambda)**I;
C     R*4   ERROR(15)  ERROR matrix (taken from TRK(I3+20) et seq.);
C     R*4   RM          = sqrt(xm**2 + ym**2);
C     I*4   MFIT        = Least-significant byte of MFITWD;
C     I*4   MZFIT       = Next-to-least-significant byte of MFITWD;
C     I*4   IFIT(2)     R-phi layer maps = ITRK(I3+35), ITRK(I3+36);
CLHB        IFIT(2)     ITRK(I3+37), ITRK(I3+38)FOR BES
C     I*4   CD          CD CODE = ITRK(I3+40).
C
C
C  Modifications:
C  ---------------------------------------------------------------------------
C  ??????????? Lan Huibin    Modified for use in BES
C
C  Jan 13,1994 Bruce Lowery  Moved WAITER, WAITERG, WAITER0, and DCRHOL
C                            into separate files.
C
C   FIX ONE BUG REALTED TO MFIT=-19 AND -9,POINTED OUT BY LIXN AND B.LOWERY
C                                              MAY-26,1994
C------------------------------------------------------------------------------