Module helixfit_track_mod

module helixfit_track_mod

        ! Uses
    use precision_mod
    use namelist_mod
    use physicalconstants_mod
    use track_mod
    use helixfit_common_mod
    use helixfit_fitpar_mod
    use helixfit_hits_mod
    use helixfit_cell_mod

        ! Types
    public type residuals_type
    public type chisqr_type
    public type kinks_type

        ! Variables
    integer, public, PARAMETER :: kAmbigModeDefault = 0
    integer, public, PARAMETER :: kAmbigModeMC = 1
    integer, public, PARAMETER :: kAmbigModeWant = 2
    integer, public, PARAMETER :: kAmbigModeSet = 3
    integer, public, PARAMETER :: kAmbigModeWC = 4
    integer, public, PARAMETER :: kResTypeNone = 0
    integer, public, PARAMETER :: kResTypeWC = 1
    integer, public, PARAMETER :: kResTypeDrift = 2
    integer, public, PARAMETER :: kResTypeTime = 3
    integer, public, PARAMETER :: kKinkTypeNone = 0
    integer, public, PARAMETER :: kKinkTypeFirst = 1
    integer, public, PARAMETER :: kKinkTypeMiddle = 2
    integer, public, PARAMETER :: kKinkTypeLast = 3
    integer, public, PARAMETER :: kKinkTypeTarget = 4
    integer, public, PARAMETER :: kKinkTypeTargetLoss = 5
    integer, public, PARAMETER :: kMaxResiduals = 200
    real (kind=R8), private, EXTERNAL :: conlev
    integer, private, EXTERNAL :: inewunit
    integer, private, EXTERNAL :: solve_linear

        ! Subroutines and functions
    public subroutine makeChiSqr (hits, res, chisqr)
    public function cmpChiSqr (a, b) result (icmp)
    private function setPlanePositions (hits, res, iplane, iresmap, V1, P1, V2, P2, V3, P3, Q, K, enableSTR, tref, tof) result (ierror)
    private subroutine goThroughFoil (istepdir, lossFoil, v, p, q)
    private subroutine setPlaneAmbig (hits, iplane, ires, res)
    private subroutine setPlaneWindows1 (hits, res, iplane, iresmap, isort, enterPos, midPos, exitPos)
    private subroutine FindClosestCluster (hits, res, iplane, iresmap, pos, firstwire, lastwire)
    private subroutine setPlaneWindows2 (hits, res, iplane, iresmap, isort, enterPos, enterMom, midPos, exitPos, exitMom)
    public function setPositions (hits, kinks, lastplane, fitpar, enableSTR, res) result (ierror)
    private subroutine doKink (kinktype, ikinkpar, fitpar, v, p, q)
    private subroutine setMcData (res)
    private subroutine sortplanehits (hits, iplane, v, isort)
    private subroutine addResidual (hits, iplane, ihit, v, p, cellpos, d, tmin, tmean, trms, isInside, wirepos, tref, tof, res, ires)
    private subroutine addKinkResidual (v, p, q, mass, kink, kinkxx0, res, ires, kinktype)
    private subroutine FindCluster (hits, iplane, iwire, firstwire, lastwire)
    private subroutine setResidualsWC (hits, fitpar, res)
    public subroutine setAmbig (hits, fitpar, res, cutoff)
    public subroutine setResidualsDrift (hits, fitpar, res)
    private subroutine setResidualsTime (hits, fitpar, res)
    private subroutine RemoveOutliers (res)
    public function removeOutliers2 (hits, res, chi) result (numRemoved)
    private function removeOutliers3 (hits, res, chi, threshold) result (numRemoved)
    public subroutine setResiduals (hits, fitpar, res, type)
    public subroutine makeTofCorrection (hits, cost, tof)
    private function findHitMC (iplane, iwire) result (ifound)
    private subroutine setDriftsMC (hits, res, tref)
    public subroutine DisableRes (res, ires)
    public subroutine ClearRes (res)
    public subroutine PrintRes (res)
    public subroutine CopyRes (dest, src)
    private function CountRes (res) result (count)
    public subroutine assignResiduals (hits, kinks)
    private function finite (name, x)
    private function V3mag (a)
    private function sqr (a) result (b)
    private function myReal8 (a4) result (a8)
    private function jsign (a) result (ival)

end module helixfit_track_mod

Description of Types

public type residuals_type
    logical :: usesSTR
             if positions use the STRs
    integer :: resType
             residuals type
    integer :: numpar
             number of fit parameters
    real (kind=R8), dimension (kMaxPar) :: par
             fit parameters (1:numpar)
    integer :: numRes
             number of residuals
    integer :: firstPlane
             first plane used
    integer :: lastPlane
             last plane used
    real (kind=R8) :: kinkfudge
             set by "setPositions()"
    integer, dimension (kMaxResiduals) :: iseq
             set by "setPositions->addResidual"
 residual sequencial number
    real (kind=R8), dimension (kMaxResiduals,2) :: trkcellpos
             track cell position
    real (kind=R8), dimension (kMaxResiduals,3) :: trkpos
             track position
    real (kind=R8), dimension (kMaxResiduals,3) :: trkmom
             track momentum
    real (kind=R8), dimension (kMaxResiduals,3) :: wirepos
             hit wire position
    integer, dimension (kMaxResiduals) :: fgindex
             FGresult%dchit() hit index
    integer, dimension (kMaxResiduals) :: dctdcindex
    integer, dimension (kMaxResiduals) :: iplane
             plane number and
    integer, dimension (kMaxResiduals) :: ihit
             hits%chits(iplane,ihit) index
    integer, dimension (kMaxResiduals) :: iwire
             wire number
    real (kind=R8), dimension (kMaxResiduals) :: dtrk
             space track
    logical, dimension (kMaxResiduals) :: isInside
             hit is inside the drift cell
    real (kind=R8), dimension (kMaxResiduals) :: ttrkmin
             set by "setPositions->addResidual"
 time track, minimum
    real (kind=R8), dimension (kMaxResiduals) :: ttrkmean
             time track, mean
    real (kind=R8), dimension (kMaxResiduals) :: ttrkrms
             time track, rms
    real (kind=R8), dimension (kMaxResiduals) :: ttrk
             time track, as used in fit
    real (kind=R8), dimension (kMaxResiduals) :: dhit
             space hit
    real (kind=R8), dimension (kMaxResiduals) :: thit
             time hit
    real (kind=R8), dimension (kMaxResiduals) :: dres
             space residual, cm
    real (kind=R8), dimension (kMaxResiduals) :: tres
             time residual, ns
    integer, dimension (kMaxResiduals) :: wantambig
             desired ambiguity: -1, +1
    integer, dimension (kMaxResiduals) :: clusterwires
             set by "setPositions->setPlanePositions->setPlaneWindows"
    real (kind=R8), dimension (kMaxResiduals) :: clustercenter
    real (kind=R8), dimension (kMaxResiduals) :: clusterwindow
    integer, dimension (kMaxResiduals) :: mcambig
             set by "setPositions->setMCdata"
 MCSP ambiguity: -1, +1
    real (kind=R8), dimension (kMaxResiduals) :: mcthit
             MCSP drift time
    real (kind=R8), dimension (kMaxResiduals) :: mcdhit
             MCSP drift distance
    real (kind=R8), dimension (kMaxResiduals) :: mcpmag
             MCSP momentum
    real (kind=R8), dimension (kMaxResiduals) :: mccost
             MCSP cos(theta)
    integer, dimension (kMaxResiduals) :: setambig
             set by "setAmbiguities"
 resolved ambiguity: -1, +1
    integer, dimension (kMaxResiduals) :: wcambig
             set by "iterate_llsq"
 wire-centers ambiguity: -1, +1
    real (kind=R8), dimension (kMaxResiduals) :: mtrk
             set by "setResiduals" according to kResTypeXXX
    real (kind=R8), dimension (kMaxResiduals) :: mhit
    real (kind=R8), dimension (kMaxResiduals) :: merr
    real (kind=R8), dimension (kMaxResiduals,3) :: hitpos
             hit position
    real (kind=R8), dimension (kMaxResiduals) :: dofIncr
             degrees of freedom increment
    real (kind=R8), dimension (kMaxResiduals) :: weights
    real (kind=R8), dimension (kMaxResiduals) :: chiSqrIncr
    real (kind=R8), dimension (kMaxResiduals) :: confLevel
    logical, dimension (kMaxResiduals) :: used
    integer, dimension (kMaxResiduals) :: ambig
    real (kind=R8), dimension (kMaxResiduals) :: clusterLikelihood
end type residuals_type

public type chisqr_type
    integer :: numhits
             number of hits
    integer :: numdof
             degrees of freedom
    real (kind=R8) :: chiSqr
             chi squared sum
    real (kind=R8) :: conflevel
             confidence level
    real (kind=R8) :: chebyshevNorm
             Chebyshev norm using narrow windows
end type chisqr_type

public type kinks_type
    real (kind=R8), dimension (kMaxPlanes) :: z
             kink z coordinate after given plane
    real (kind=R8), dimension (kMaxPlanes) :: radlen
             kink radiation length per unit
    integer, dimension (kMaxPlanes) :: kinktype
             kink type: see kKinkTypeXXX
    integer, dimension (kMaxPlanes,kMaxHitsPerPlane) :: ihitres
             hit residual index into residuals_type
    integer, dimension (kMaxPlanes) :: ikinkres
             kink residual index into residuals_type
    integer, dimension (kMaxPlanes) :: ikinkpar
             index to the kink parameters in fitpar%par
end type kinks_type

Description of Variables

integer, public, PARAMETER :: kAmbigModeDefault = 0
 HelixFitAmbigMode values:

integer, public, PARAMETER :: kAmbigModeMC = 1
 take ambiguities from MCSP data

integer, public, PARAMETER :: kAmbigModeWant = 2
 set ambiguities as they wish

integer, public, PARAMETER :: kAmbigModeSet = 3
 use resolved ambiguities

integer, public, PARAMETER :: kAmbigModeWC = 4
 use wire-centers ambiguities

integer, public, PARAMETER :: kResTypeNone = 0
 before call to setResiduals()

integer, public, PARAMETER :: kResTypeWC = 1
 wire-center residuals

integer, public, PARAMETER :: kResTypeDrift = 2
 drift distance residuals

integer, public, PARAMETER :: kResTypeTime = 3
 drift time residuals

integer, public, PARAMETER :: kKinkTypeNone = 0
 No kink

integer, public, PARAMETER :: kKinkTypeFirst = 1
 First kink

integer, public, PARAMETER :: kKinkTypeMiddle = 2
 normal kink

integer, public, PARAMETER :: kKinkTypeLast = 3
 Last kink

integer, public, PARAMETER :: kKinkTypeTarget = 4
 Target kink

integer, public, PARAMETER :: kKinkTypeTargetLoss = 5
 Target kink with energy loss

integer, public, PARAMETER :: kMaxResiduals = 200
 maximum number of hits

real (kind=R8), private, EXTERNAL :: conlev
 external confidence level calculator

integer, private, EXTERNAL :: inewunit
 external i/o unit allocator

integer, private, EXTERNAL :: solve_linear
 external solve linear equations

Description of Subroutines and Functions

public subroutine makeChiSqr (hits, res, chisqr)
    type (hits_type), INTENT(in) :: hits
            
 Description: compute ChiSqr and confidence level

 Arguments:

 hits
    type (residuals_type), INTENT(in) :: res
             input residuals
    type (chisqr_type), INTENT(out) :: chisqr
             output chi squared
    ! Calls: PrintRes, die
end subroutine makeChiSqr

public function cmpChiSqr (a, b) result (icmp)
    type (chisqr_type), INTENT(in) :: a
            
 Description: compare two confidence levels:
    the larger conf.level (smaller chiSqr/numdof)
    is the better.
 Return:
    -1: "a" is worse than "b"
     0: "a" and "b" are equal
     1: "a" is better than "b"
    type (chisqr_type), INTENT(in) :: b
    integer :: icmp
             return value
end function cmpChiSqr

private function setPlanePositions (hits, res, iplane, iresmap, V1, P1, V2, P2, V3, P3, Q, K, enableSTR, tref, tof) result (ierror)
    type (hits_type), INTENT(in) :: hits
            
 Description: calculate track positions along the trajectory

 Arguments:

 hits
    type (residuals_type), INTENT(inout) :: res
             output residuals
    integer, INTENT(in) :: iplane
    integer, INTENT(in), dimension (:) :: iresmap
    real (kind=R8), INTENT(in), dimension (3) :: V1
             at plane entrance
    real (kind=R8), INTENT(in), dimension (3) :: P1
             at plane entrance
    real (kind=R8), INTENT(in), dimension (3) :: V2
             at plane z
    real (kind=R8), INTENT(in), dimension (3) :: P2
             at plane z
    real (kind=R8), INTENT(in), dimension (3) :: V3
             at plane exit
    real (kind=R8), INTENT(in), dimension (3) :: P3
             at plane exit
    integer, INTENT(in) :: Q
    real (kind=R8), INTENT(in) :: K
    logical, INTENT(in) :: enableSTR
    real (kind=R8), INTENT(in) :: tref
    real (kind=R8), INTENT(in) :: tof
             time of flight correction
    integer :: ierror
             return value
    ! Calls: addresidual, setPlaneWindows2, sortplanehits
end function setPlanePositions

private subroutine goThroughFoil (istepdir, lossFoil, v, p, q)
    integer, INTENT(in) :: istepdir
            
 Lose energy when going through a plane foil

 stepping direction: -1 or +1
    real (kind=R8), INTENT(in) :: lossFoil
    real (kind=R8), INTENT(inout), dimension (3) :: v
    real (kind=R8), INTENT(inout), dimension (3) :: p
    integer, INTENT(inout) :: q
end subroutine goThroughFoil

private subroutine setPlaneAmbig (hits, iplane, ires, res)
    type (hits_type), INTENT(in) :: hits
            
 Description: resolve drift ambiguities using plane hit multiplicity

 Arguments:
    integer, INTENT(in) :: iplane
    integer, INTENT(in), dimension (:) :: ires
    type (residuals_type), INTENT(inout) :: res
    ! Calls: setHitAmbig
end subroutine setPlaneAmbig

private subroutine setPlaneWindows1 (hits, res, iplane, iresmap, isort, enterPos, midPos, exitPos)
    type (hits_type), INTENT(in) :: hits
            
 Description: set residuals cluster positions

 Arguments:
    type (residuals_type), INTENT(inout) :: res
    integer, INTENT(in) :: iplane
    integer, INTENT(in), dimension (:) :: iresmap
    integer, INTENT(in), dimension (:) :: isort
    real (kind=R8), INTENT(in), dimension (3) :: enterPos
             coordinates at enter foil plane
    real (kind=R8), INTENT(in), dimension (3) :: midPos
             coordinates at wires plane
    real (kind=R8), INTENT(in), dimension (3) :: exitPos
             coordinates at exit foil plane
    ! Calls: FindCluster
end subroutine setPlaneWindows1

private subroutine FindClosestCluster (hits, res, iplane, iresmap, pos, firstwire, lastwire)
    type (hits_type), INTENT(in) :: hits
    type (residuals_type), INTENT(in) :: res
    integer, INTENT(in) :: iplane
    integer, INTENT(in), dimension (:) :: iresmap
    real (kind=R8), INTENT(in), dimension (3) :: pos
    integer, INTENT(out) :: firstwire
    integer, INTENT(out) :: lastwire
    ! Calls: FindCluster
end subroutine FindClosestCluster

private subroutine setPlaneWindows2 (hits, res, iplane, iresmap, isort, enterPos, enterMom, midPos, exitPos, exitMom)
    type (hits_type), INTENT(in) :: hits
            
 Description: set residuals cluster positions

 Arguments:
    type (residuals_type), INTENT(inout) :: res
    integer, INTENT(in) :: iplane
    integer, INTENT(in), dimension (:) :: iresmap
    integer, INTENT(in), dimension (:) :: isort
    real (kind=R8), INTENT(in), dimension (3) :: enterPos
             coordinates at enter foil plane
    real (kind=R8), INTENT(in), dimension (3) :: enterMom
             momentum at enter foil plane
    real (kind=R8), INTENT(in), dimension (3) :: midPos
             coordinates at wires plane
    real (kind=R8), INTENT(in), dimension (3) :: exitPos
             coordinates at exit foil plane
    real (kind=R8), INTENT(in), dimension (3) :: exitMom
             momentum at exit foil plane
    ! Calls: FindClosestCluster
end subroutine setPlaneWindows2

public function setPositions (hits, kinks, lastplane, fitpar, enableSTR, res) result (ierror)
    type (hits_type), INTENT(in) :: hits
            
 Description: calculate track positions along the trajectory

 Arguments:

 hits
    type (kinks_type), INTENT(in) :: kinks
             kinks
    integer, INTENT(in) :: lastplane
             do not go past this plane
    type (fitpar_type), INTENT(in) :: fitpar
             input track
    logical, INTENT(in) :: enableSTR
             enable STR calculations
    type (residuals_type), INTENT(out) :: res
             output residuals
    integer :: ierror
            
 return value
    ! Calls: ClearRes, addKinkResidual, die, doKink, goThroughFoil, kerror, makeTofCorrection, setMcData, setPlaneAmbig
end function setPositions

private subroutine doKink (kinktype, ikinkpar, fitpar, v, p, q)
    integer, INTENT(in) :: kinktype
            
 Apply the kink to the track position and momentum

 kink type kKinkTypeXXX

Author: R. Bayes
    integer, INTENT(in) :: ikinkpar
             offset to kink parameters
    type (fitpar_type), INTENT(in) :: fitpar
             fit parameters
    real (kind=R8), INTENT(inout), dimension (3) :: v
             track position
    real (kind=R8), INTENT(inout), dimension (3) :: p
             track momentum
    integer, INTENT(inout) :: q
             track charge
end subroutine doKink

private subroutine setMcData (res)
    type (residuals_type), INTENT(inout) :: res
            
 Description: set the MC ambiguity from UnpMC_mod::MCSP

 Arguments:

 residuals
end subroutine setMcData

private subroutine sortplanehits (hits, iplane, v, isort)
    type (hits_type), INTENT(in) :: hits
            
 Description: Sort hits on the plane in order of increasing
 distance from point v()

 hits
    integer, INTENT(in) :: iplane
    real (kind=R8), INTENT(in), dimension (3) :: v
    integer, INTENT(out), dimension (kMaxHitsPerPlane) :: isort
end subroutine sortplanehits

private subroutine addResidual (hits, iplane, ihit, v, p, cellpos, d, tmin, tmean, trms, isInside, wirepos, tref, tof, res, ires)
    type (hits_type), INTENT(in) :: hits
             hits record
    integer, INTENT(in) :: iplane
             hit index in
    integer, INTENT(in) :: ihit
             hits%chits(iplane,ihit)
    real (kind=R8), INTENT(in), dimension (3) :: v
             track position
    real (kind=R8), INTENT(in), dimension (3) :: p
             track momentum
    real (kind=R8), INTENT(in), dimension (2) :: cellpos
             track cell position
    real (kind=R8), INTENT(in) :: d
             track drift distance
    real (kind=R8), INTENT(in) :: tmin
             track drift time
    real (kind=R8), INTENT(in) :: tmean
             track drift time, clustered mean
    real (kind=R8), INTENT(in) :: trms
             track drift time, clustered rms
    logical, INTENT(in) :: isInside
             hit is inside the drift cell
    real (kind=R8), INTENT(in), dimension (3) :: wirepos
             hit wire position
    real (kind=R8), INTENT(in) :: tref
             track start time
    real (kind=R8), INTENT(in) :: tof
             time of flight correction
    type (residuals_type), INTENT(inout) :: res
             residuals
    integer, INTENT(in) :: ires
             residual index
    ! Calls: DisableRes
end subroutine addResidual

private subroutine addKinkResidual (v, p, q, mass, kink, kinkxx0, res, ires, kinktype)
    real (kind=R8), INTENT(in), dimension (3) :: v
             track position
    real (kind=R8), INTENT(in), dimension (3) :: p
             track momentum
    integer, INTENT(in) :: q
             track charge
    real (kind=R8), INTENT(in) :: mass
             particle mass in MeV/c^2
    real (kind=R8), INTENT(in) :: kink
             track kink
    real (kind=R8), INTENT(in) :: kinkxx0
             kink radiation length x/x0
    type (residuals_type), INTENT(inout) :: res
             residuals
    integer, INTENT(in) :: ires
             residual index
    integer, INTENT(in) :: kinktype
             kink type
end subroutine addKinkResidual

private subroutine FindCluster (hits, iplane, iwire, firstwire, lastwire)
    type (hits_type), INTENT(in) :: hits
            
 Description: return first and last wire of the hits cluster
              that includes "iwire" in plane "iplane", 

 hits record
    integer, INTENT(in) :: iplane
    integer, INTENT(in) :: iwire
    integer, INTENT(out) :: firstwire
    integer, INTENT(out) :: lastwire
    ! Calls: die
end subroutine FindCluster

private subroutine setResidualsWC (hits, fitpar, res)
    type (hits_type), INTENT(in) :: hits
             hits record
    type (fitpar_type), INTENT(in) :: fitpar
    type (residuals_type), INTENT(inout) :: res
             residuals
end subroutine setResidualsWC

public subroutine setAmbig (hits, fitpar, res, cutoff)
    type (hits_type), INTENT(in) :: hits
            
 Description: Set ambiguities

 hits record
    type (fitpar_type), INTENT(in) :: fitpar
             fit parameters
    type (residuals_type), INTENT(inout) :: res
             residuals
    real (kind=R8), INTENT(in) :: cutoff
    ! Calls: die
end subroutine setAmbig

public subroutine setResidualsDrift (hits, fitpar, res)
    type (hits_type), INTENT(in) :: hits
             hits record
    type (fitpar_type), INTENT(in) :: fitpar
             fit parameters
    type (residuals_type), INTENT(inout) :: res
             residuals
    ! Calls: disableRes
end subroutine setResidualsDrift

private subroutine setResidualsTime (hits, fitpar, res)
    type (hits_type), INTENT(in) :: hits
    type (fitpar_type), INTENT(in) :: fitpar
    type (residuals_type), INTENT(inout) :: res
end subroutine setResidualsTime

private subroutine RemoveOutliers (res)
    type (residuals_type), INTENT(inout) :: res
    ! Calls: sortreal8
end subroutine RemoveOutliers

public function removeOutliers2 (hits, res, chi) result (numRemoved)
    type (hits_type), INTENT(in) :: hits
            
 Driver for outlier removal. Calls removeOutliers3()
    type (residuals_type), INTENT(inout) :: res
    type (chisqr_type), INTENT(inout) :: chi
    integer :: numRemoved
             return value
end function removeOutliers2

private function removeOutliers3 (hits, res, chi, threshold) result (numRemoved)
    type (hits_type), INTENT(in) :: hits
            
 Remove hits with confLevel below given threshold
    type (residuals_type), INTENT(inout) :: res
    type (chisqr_type), INTENT(inout) :: chi
    real (kind=R8), INTENT(in) :: threshold
             removal threshold
    integer :: numRemoved
             return value
    ! Calls: makeChiSqr
end function removeOutliers3

public subroutine setResiduals (hits, fitpar, res, type)
    type (hits_type), INTENT(in) :: hits
    type (fitpar_type), INTENT(in) :: fitpar
    type (residuals_type), INTENT(inout) :: res
    integer, INTENT(in) :: type
             kResTypeXXX
    ! Calls: RemoveOutliers, die, setResidualsDrift, setResidualsTime, setResidualsWC
end subroutine setResiduals

public subroutine makeTofCorrection (hits, cost, tof)
    type (hits_type), INTENT(in) :: hits
            
 Description: set the time-of-flight correction

 hits
    real (kind=R8), INTENT(in) :: cost
             cos(theta)
    real (kind=R8), INTENT(out), dimension (kMaxPlanes) :: tof
             per-plane tof correction
end subroutine makeTofCorrection

private function findHitMC (iplane, iwire) result (ifound)
    integer, INTENT(in) :: iplane
    integer, INTENT(in) :: iwire
    integer :: ifound
             return value
end function findHitMC

private subroutine setDriftsMC (hits, res, tref)
    type (hits_type), INTENT(in) :: hits
    type (residuals_type), INTENT(inout) :: res
    real (kind=R8), INTENT(in) :: tref
    ! Calls: DisableRes
end subroutine setDriftsMC

public subroutine DisableRes (res, ires)
    type (residuals_type), INTENT(inout) :: res
    integer, INTENT(in) :: ires
end subroutine DisableRes

public subroutine ClearRes (res)
    type (residuals_type), INTENT(out) :: res
end subroutine ClearRes

public subroutine PrintRes (res)
    type (residuals_type), INTENT(in) :: res
end subroutine PrintRes

public subroutine CopyRes (dest, src)
    type (residuals_type), INTENT(out) :: dest
    type (residuals_type), INTENT(in) :: src
end subroutine CopyRes

private function CountRes (res) result (count)
    type (residuals_type), INTENT(in) :: res
            
 Description: return the number of "used" residuals
    integer :: count
             return value
end function CountRes

public subroutine assignResiduals (hits, kinks)
    type (hits_type), INTENT(in) :: hits
    type (kinks_type), INTENT(inout) :: kinks
end subroutine assignResiduals

private function finite (name, x)
    character (len=*), INTENT(in) :: name
            
 Description: return .true. if (x) is not a NaN
    real (kind=R8), INTENT(in) :: x
    logical :: finite
    ! Calls: kerror
end function finite

private function V3mag (a)
    real (kind=R8), INTENT(in), dimension (3) :: a
            
 Description: compute vector length
    real (kind=R8) :: v3mag
end function V3mag

private function sqr (a) result (b)
    real (kind=R8), INTENT(in) :: a
             Description: return square of a real(r8) number
    real (kind=R8) :: b
end function sqr

private function myReal8 (a4) result (a8)
    real, INTENT(in) :: a4
    real (kind=R8) :: a8
end function myReal8

private function jsign (a) result (ival)
    real (kind=R8), INTENT(in) :: a
            
 Description: return sign of a real(r8) number
 Returns: +1 if positive
          -1 if negative
           0 if zero
    integer :: ival
end function jsign