Module firstguessko_mod

module firstguessko_mod

        ! Uses
    use precision_mod
    use namelist_mod
    use physicalconstants_mod, only: pi
    use det_geom_mod, ONLY: max_planes_d, max_planes_p
    use timezero_mod, ONLY:timezero_costheta

        ! Types
    private type hits_type
    private type uvfit_type
    private type zphifit_type
    private type fit_type

        ! Variables
    integer, private, PARAMETER :: kMaxHits = max_planes_d / 2 + max_planes_p / 2
    integer, private, PARAMETER :: kMaxMult = 10
    integer, private, PARAMETER :: kMaxTdc = 100
    integer, private, PARAMETER :: kMinHits = 4
    real, private, PARAMETER :: kMinZ = 1.0
    real, private, PARAMETER :: kMaxZ = 20.0
    real, private, PARAMETER :: kMinCost = 0.01
    real, private, PARAMETER :: kMinRadius = 0.3000
    real, private, PARAMETER :: kMinMuonRadius = 0.05
    real, private, PARAMETER :: kMaxRadius = 15.0
    real, private, PARAMETER :: kMaxTant = 7.999
    real, private, PARAMETER :: kCutQ = 10.0
    real, private, PARAMETER :: kCutRQ = 0.4
    real, private, PARAMETER :: kHUGEVAL = 1E37
    real, private, PARAMETER :: kRadiusWindingGuessChanged = 1.
    integer, private, PARAMETER :: kSuccess = 0
    integer, private, PARAMETER :: kFailUnknown = 1
    integer, private, PARAMETER :: kFailMinHits = 2
    integer, private, PARAMETER :: kFailTooShort = 3
    integer, private, PARAMETER :: kFailZeroRadius = 4
    integer, private, PARAMETER :: kFailMaxTant = 5
    integer, private, PARAMETER :: kFailBadRadius = 6
    integer, private, PARAMETER :: kFailAnchors = 7
    integer, private, PARAMETER :: kFailBadPhaseCost = 8
    integer, private, PARAMETER :: kFailBadPhaseTant = 9
    integer, private, PARAMETER :: kFailNoCircles = 10
    integer, private, PARAMETER :: kFailNoTrack = 11
    integer, private, PARAMETER :: kFailHasGap = 12
    integer (kind=I4), private, EXTERNAL :: iutime

        ! Subroutines and functions
    public subroutine FirstGuessKOwindow (iwindow, iStream)
    private subroutine doGeneric (iwindow, iStream)
    private subroutine removeHits (hits, fit)
    private subroutine doMuon (iwindow)
    public subroutine FirstGuessKOmuon (iwindow, fgresult)
    private subroutine NewUVfit (uvfit)
    private subroutine NewZPhiFit (zphifit)
    private subroutine NewFit (fit)
    private subroutine PrintZPhiFit (hits, fit, doPrintHits)
    private subroutine storeFit (hits, fit, muonmode, fg)
    private subroutine addFit (hits, fit, muonmode)
    private subroutine findHits (iwindow, iStream, muonmode, hits)
    private subroutine ClearHits (hits)
    private subroutine storeHit (hits, isDC, isFullyInstrumented, u, v, du, dv, z, s, w, numtdchits, tdchits)
    private subroutine MergeSort (hits, hits1, hits2)
    private subroutine copyHit (dst, idst, src, isrc)
    private subroutine printHits (hits)
    private subroutine iterateSubtracks (hits, muonmode, fit)
    private subroutine iterateUV (hits, ifirstHit, ilastHit, fit, muonmode)
    private subroutine FudgeDeltaTrack (uv, phase)
    private subroutine ComputeEdges (hits, ifirsthit, ilasthit, isEdge)
    private subroutine ComputeAverageMultiplicity (hits, iFirstHit, iLastHit, avgU, avgV, avgA)
    private subroutine ComputeTant (hits, iFirstHit, iLastHit, minTant, maxTant)
    private function FindMuonCircles (hits, iFirstHit, iLastHit, idir, isEdge, numuvc, uvc, uvcindex)
    private function FindCircles (hits, ifirsthit, ilasthit, idir, isEdge, numuvc, uvc, uvcindex)
    private subroutine FitUV5 (hits, ia1, im1, ia2, im2, ia3, im3, fit)
    private subroutine makeCircle (u, v, uc, vc, rc)
    private subroutine ComputeCircleNorm (hits, ifirst, ilast, isEdge, fituv)
    private subroutine iteratePhase (hits, iFirstHit, iLastHit, isEdge, fit, uv, muonmode, phase)
    private subroutine tryPhase (hits, hitsr, hitsp, hitsdp, r0, ifirst, ilast, ia1, im1, ia2, im2, iwind, minTant, maxTant, muonmode, phase)
    private function angleDiff (a, b)
    private function reduce (a)
    private function rad2deg (a)
    private function deg2rad (a)
    private function jsign (a) result (ival)
    private function sqr (x) result (x2)
    private subroutine Fgkodie ()
    private function kcutrq2 (x) result (x2)

end module firstguessko_mod

Description of Types

private type hits_type
    integer :: iwindow
             this window
    integer :: numHits
             total number of hits
    integer, dimension (kMaxHits) :: mult
             number of clusters per hit
    logical, dimension (kMaxHits) :: isDChit
             TRUE for DC hits, FALSE for others
    logical, dimension (kMaxHits) :: isFullyInstrumented
             TRUE for fully instrumented planes
    real (kind=R8), dimension (kMaxHits) :: z
             cluster z, cm
    real (kind=R8), dimension (kMaxHits,kMaxMult) :: u
             cluster center, cm
    real (kind=R8), dimension (kMaxHits,kMaxMult) :: v
             cluster center, cm
    real (kind=R8), dimension (kMaxHits,kMaxMult) :: du
             cluster u size, cm
    real (kind=R8), dimension (kMaxHits,kMaxMult) :: dv
             cluster v size, cm
    real (kind=R8), dimension (kMaxHits,kMaxMult) :: s
             cluster uv diagonal, cm
    real (kind=R8), dimension (kMaxHits,kMaxMult) :: w
             cluster weight?
    integer, dimension (kMaxHits,kMaxMult,0:kMaxTdc) :: tdchits
             indices of TDC hits
end type hits_type

private type uvfit_type
    real (kind=R8) :: u
             circle center and radius
    real (kind=R8) :: v
             circle center and radius
    real (kind=R8) :: r
             circle center and radius
    integer :: numhits
             number of hits used in fit
    logical, dimension (kMaxHits) :: isGood
             save only good hits
    integer, dimension (kMaxHits) :: ihit
             index into hits_type%z(), etc
    integer, dimension (kMaxHits) :: imult
             index into hits_type%u(ihit,imult), etc
    real (kind=R8), dimension (kMaxHits,kMaxMult) :: nres
             norm for each hit
    integer :: ia1
    integer :: ia2
    integer :: ia3
    integer :: im1
    integer :: im2
    integer :: im3
    real (kind=R8) :: norm
             Chebyshev norm
end type uvfit_type

private type zphifit_type
    integer :: ifgerror
             error code
    integer :: ifirst
             first and last hits of subtrack
    integer :: ilast
             first and last hits of subtrack
    integer :: ianchor1
             helix anchors
    integer :: imult1
             helix anchors
    integer :: ianchor2
             helix anchors
    integer :: imult2
             helix anchors
    integer :: iwind
             winding number
    real (kind=R8) :: avgwind
    real (kind=R8) :: avgcost
    real (kind=R8) :: phi0
             phase at z=0
    real (kind=R8) :: afreq
             angular frequency
    real (kind=R8) :: tant
             tan(theta)
    real (kind=R8) :: cost
             cos(theta)
    integer :: numhits
             number of hits used in fit
    integer, dimension (kMaxHits) :: ihit
             index into hits_type%z(), etc
    integer, dimension (kMaxHits) :: imult
             index into hits_type%u(ihit,imult), etc
    logical, dimension (kMaxHits) :: isGood
             save only good hits
    real (kind=R8), dimension (kMaxHits,kMaxMult) :: pwidth
             width of each hit
    real (kind=R8), dimension (kMaxHits,kMaxMult) :: pres
             angle residual for each hit
    real (kind=R8), dimension (kMaxHits,kMaxMult) :: qres
             quality function for each hit
    real (kind=R8), dimension (kMaxHits,kMaxMult) :: rres
             radius residual for each hit
    real (kind=R8) :: qfunc
             quality function
end type zphifit_type

private type fit_type
    integer :: ifirsthit
             first and last hits used by this fit
    integer :: ilasthit
             first and last hits used by this fit
    integer :: idir
             direction we should loop over hits
    integer :: ifgerror
             error code for this fit
    type (uvfit_type) :: uv
             u-v circle fit for this fit
    real (kind=R8) :: clusterSize
             average cluster size
    real (kind=R8) :: minTant
             minimum tan(theta) from cluster size 
    real (kind=R8) :: maxTant
             maximum tan(theta) from cluster size 
    type (zphifit_type) :: zphi
             phase fit
end type fit_type

Description of Variables

integer, private, PARAMETER :: kMaxHits = max_planes_d / 2 + max_planes_p / 2

integer, private, PARAMETER :: kMaxMult = 10
 maximum cluster multiplicity

integer, private, PARAMETER :: kMaxTdc = 100
 maximum number of hits per cluster

integer, private, PARAMETER :: kMinHits = 4
 min number of hits to make a track

real, private, PARAMETER :: kMinZ = 1.0

real, private, PARAMETER :: kMaxZ = 20.0
 maximum distance between hits, cm

real, private, PARAMETER :: kMinCost = 0.01
 minimum cos(theta)

real, private, PARAMETER :: kMinRadius = 0.3000
 minimum radius, cm

real, private, PARAMETER :: kMinMuonRadius = 0.05
 min radius for muonmode, cm

real, private, PARAMETER :: kMaxRadius = 15.0
 maximum radius, cm

real, private, PARAMETER :: kMaxTant = 7.999
 maximum tan(theta)

real, private, PARAMETER :: kCutQ = 10.0
 degrees

real, private, PARAMETER :: kCutRQ = 0.4
 cm

real, private, PARAMETER :: kHUGEVAL = 1E37
 huge value

real, private, PARAMETER :: kRadiusWindingGuessChanged = 1.
cm

integer, private, PARAMETER :: kSuccess = 0
 no error

integer, private, PARAMETER :: kFailUnknown = 1
 unknown, unexpected error

integer, private, PARAMETER :: kFailMinHits = 2
 too few hits

integer, private, PARAMETER :: kFailTooShort = 3
 track is too short in z

integer, private, PARAMETER :: kFailZeroRadius = 4
 helix radius is zero

integer, private, PARAMETER :: kFailMaxTant = 5
 failed max Tan(theta) cut

integer, private, PARAMETER :: kFailBadRadius = 6
 helix radius is too big

integer, private, PARAMETER :: kFailAnchors = 7
 cannot find helix anchors

integer, private, PARAMETER :: kFailBadPhaseCost = 8
 xxx

integer, private, PARAMETER :: kFailBadPhaseTant = 9
 xxx

integer, private, PARAMETER :: kFailNoCircles = 10
 cannot find any circle candidates

integer, private, PARAMETER :: kFailNoTrack = 11
 cannot find any tracks

integer, private, PARAMETER :: kFailHasGap = 12
 track has a gap

integer (kind=I4), private, EXTERNAL :: iutime

Description of Subroutines and Functions

public subroutine FirstGuessKOwindow (iwindow, iStream)
    integer, INTENT(in) :: iwindow
    integer, INTENT(IN) :: iStream
    ! Calls: doGeneric, doMuon
end subroutine FirstGuessKOwindow

private subroutine doGeneric (iwindow, iStream)
    integer, INTENT(in) :: iwindow
    integer, INTENT(IN) :: iStream
    ! Calls: addFit, findHits, iterateSubtracks, printHits, removeHits
end subroutine doGeneric

private subroutine removeHits (hits, fit)
    type (hits_type), INTENT(inout) :: hits
    type (fit_type), INTENT(in) :: fit
    ! Calls: compressHits, removeMult
end subroutine removeHits

private subroutine doMuon (iwindow)
    integer, INTENT(in) :: iwindow
    ! Calls: FirstGuessKOmuon
end subroutine doMuon

public subroutine FirstGuessKOmuon (iwindow, fgresult)
    integer, INTENT(in) :: iwindow
    type (FGresult_type), INTENT(out) :: fgresult
    ! Calls: NewFit, findHits, iterateUV, printHits, storeFit
end subroutine FirstGuessKOmuon

private subroutine NewUVfit (uvfit)
    type (uvfit_type), INTENT(out) :: uvfit
end subroutine NewUVfit

private subroutine NewZPhiFit (zphifit)
    type (zphifit_type), INTENT(out) :: zphifit
end subroutine NewZPhiFit

private subroutine NewFit (fit)
    type (fit_type), INTENT(out) :: fit
    ! Calls: NewUVFit, NewZPhiFit
end subroutine NewFit

private subroutine PrintZPhiFit (hits, fit, doPrintHits)
    type (hits_type), INTENT(in) :: hits
    type (zphifit_type), INTENT(in) :: fit
    logical, INTENT(in) :: doPrintHits
end subroutine PrintZPhiFit

private subroutine storeFit (hits, fit, muonmode, fg)
    type (hits_type), INTENT(in) :: hits
    type (fit_type), INTENT(in) :: fit
    logical, INTENT(in) :: muonmode
    type (FGresult_type), INTENT(out) :: fg
    ! Calls: printzphifit
end subroutine storeFit

private subroutine addFit (hits, fit, muonmode)
    type (hits_type), INTENT(in) :: hits
    type (fit_type), INTENT(in) :: fit
    logical, INTENT(in) :: muonmode
    ! Calls: filltree, makebranch, setbranchvar, storeFit
end subroutine addFit

private subroutine findHits (iwindow, iStream, muonmode, hits)
    integer, INTENT(in) :: iwindow
    integer, INTENT(IN) :: iStream
    logical, INTENT(in) :: muonmode
    type (hits_type), INTENT(out) :: hits
    ! Calls: ClearHits, MergeSort, storeHit
end subroutine findHits

private subroutine ClearHits (hits)
    type (hits_type), INTENT(out) :: hits
end subroutine ClearHits

private subroutine storeHit (hits, isDC, isFullyInstrumented, u, v, du, dv, z, s, w, numtdchits, tdchits)
    type (hits_type), INTENT(inout) :: hits
    logical, INTENT(in) :: isDC
    logical, INTENT(in) :: isFullyInstrumented
    real (kind=R8), INTENT(in) :: u
    real (kind=R8), INTENT(in) :: v
    real (kind=R8), INTENT(in) :: du
    real (kind=R8), INTENT(in) :: dv
    real (kind=R8), INTENT(in) :: z
    real (kind=R8), INTENT(in) :: s
    real (kind=R8), INTENT(in) :: w
    integer, INTENT(in) :: numtdchits
    integer, INTENT(in), dimension (:) :: tdchits
end subroutine storeHit

private subroutine MergeSort (hits, hits1, hits2)
    type (hits_type), INTENT(out) :: hits
    type (hits_type), INTENT(in) :: hits1
    type (hits_type), INTENT(in) :: hits2
    ! Calls: copyHit, fgkodie
end subroutine MergeSort

private subroutine copyHit (dst, idst, src, isrc)
    type (hits_type), INTENT(inout) :: dst
    integer, INTENT(in) :: idst
    type (hits_type), INTENT(in) :: src
    integer, INTENT(in) :: isrc
end subroutine copyHit

private subroutine printHits (hits)
    type (hits_type), INTENT(in) :: hits
end subroutine printHits

private subroutine iterateSubtracks (hits, muonmode, fit)
    type (hits_type), INTENT(in) :: hits
    logical, INTENT(in) :: muonmode
    type (fit_type), INTENT(out) :: fit
    ! Calls: iterateUV
end subroutine iterateSubtracks
$    if (rootTreeType .eq. 5) then
$       if (firsttime) then
$          firsttime = .false.
$          call makebranch("FGKO", &
$               "numhits/I:"// &
$               "mean/D:rms/D:isect/I:mean1/D:mean2/D:rms1/D:rms2/D"// &
$               "")
$          call setbranchvar("FGKO","numhits",hits%numhits)
$          call setbranchvar("FGKO","mean", meanmult)
$          call setbranchvar("FGKO","rms",  rmsmult)
$          call setbranchvar("FGKO","isect",bestsect)
$          call setbranchvar("FGKO","mean1",bestmult1)
$          call setbranchvar("FGKO","mean2",bestmult2)
$          call setbranchvar("FGKO","rms1", bestrms1)
$          call setbranchvar("FGKO","rms2", bestrms2)
$       end if
$
$       call filltree()
$    end if

private subroutine iterateUV (hits, ifirstHit, ilastHit, fit, muonmode)
    type (hits_type), INTENT(in) :: hits
    integer, INTENT(in) :: iFirstHit
    integer, INTENT(in) :: iLastHit
    type (fit_type), INTENT(out) :: fit
    logical, INTENT(in) :: muonmode
    ! Calls: ComputeAverageMultiplicity, ComputeEdges, ComputeTant, FudgeDeltaTrack, NewFit, NewZPhiFit, NewZphiFit, addFit, iteratePhase
end subroutine iterateUV

private subroutine FudgeDeltaTrack (uv, phase)
    type (uvfit_type), INTENT(inout) :: uv
    type (zphifit_type), INTENT(inout) :: phase
end subroutine FudgeDeltaTrack

private subroutine ComputeEdges (hits, ifirsthit, ilasthit, isEdge)
    type (hits_type), INTENT(in) :: hits
    integer, INTENT(in) :: iFirstHit
    integer, INTENT(in) :: iLastHit
    logical, INTENT(out), dimension (kMaxHits,kMaxMult) :: isEdge
end subroutine ComputeEdges

private subroutine ComputeAverageMultiplicity (hits, iFirstHit, iLastHit, avgU, avgV, avgA)
    type (hits_type), INTENT(in) :: hits
    integer, INTENT(in) :: iFirstHit
    integer, INTENT(in) :: iLastHit
    real (kind=R8), INTENT(out) :: avgU
    real (kind=R8), INTENT(out) :: avgV
    real (kind=R8), INTENT(out) :: avgA
end subroutine ComputeAverageMultiplicity

private subroutine ComputeTant (hits, iFirstHit, iLastHit, minTant, maxTant)
    type (hits_type), INTENT(in) :: hits
    integer, INTENT(in) :: iFirstHit
    integer, INTENT(in) :: iLastHit
    real (kind=R8), INTENT(out) :: minTant
    real (kind=R8), INTENT(out) :: maxTant
end subroutine ComputeTant

private function FindMuonCircles (hits, iFirstHit, iLastHit, idir, isEdge, numuvc, uvc, uvcindex)
    type (hits_type), INTENT(in) :: hits
    integer, INTENT(in) :: iFirstHit
    integer, INTENT(in) :: iLastHit
    integer, INTENT(in) :: idir
    logical, INTENT(in), dimension (kMaxHits,kMaxMult) :: isEdge
    integer, INTENT(out) :: numuvc
    type (uvfit_type), INTENT(out), dimension (:) :: uvc
    integer, INTENT(out), dimension (:) :: uvcindex
    integer :: FindMuonCircles
             return value
    ! Calls: ComputeCircleNorm
end function FindMuonCircles

private function FindCircles (hits, ifirsthit, ilasthit, idir, isEdge, numuvc, uvc, uvcindex)
    type (hits_type), INTENT(in) :: hits
    integer, INTENT(in) :: iFirstHit
    integer, INTENT(in) :: iLastHit
    integer, INTENT(in) :: idir
    logical, INTENT(in), dimension (kMaxHits,kMaxMult) :: isEdge
    integer, INTENT(out) :: numuvc
    type (uvfit_type), INTENT(out), dimension (:) :: uvc
    integer, INTENT(out), dimension (:) :: uvcindex
    integer :: FindCircles
             return value
    ! Calls: ComputeCircleNorm, fitUV5
end function FindCircles

private subroutine FitUV5 (hits, ia1, im1, ia2, im2, ia3, im3, fit)
    type (hits_type), INTENT(in) :: hits
    integer, INTENT(in) :: ia1
    integer, INTENT(in) :: im1
    integer, INTENT(in) :: ia2
    integer, INTENT(in) :: im2
    integer, INTENT(in) :: ia3
    integer, INTENT(in) :: im3
    type (uvfit_type), INTENT(out) :: fit
    ! Calls: NewUVfit, makeCircle
end subroutine FitUV5

private subroutine makeCircle (u, v, uc, vc, rc)
    real (kind=R8), INTENT(in), dimension (3) :: u
    real (kind=R8), INTENT(in), dimension (3) :: v
    real (kind=R8), INTENT(out) :: uc
    real (kind=R8), INTENT(out) :: vc
    real (kind=R8), INTENT(out) :: rc
end subroutine makeCircle

private subroutine ComputeCircleNorm (hits, ifirst, ilast, isEdge, fituv)
    type (hits_type), intent(in) :: hits
    integer, intent(in) :: ifirst
    integer, intent(in) :: ilast
    logical, INTENT(in), dimension (kMaxHits,kMaxMult) :: isEdge
    type (uvfit_type), intent(inout) :: fituv
end subroutine ComputeCircleNorm

private subroutine iteratePhase (hits, iFirstHit, iLastHit, isEdge, fit, uv, muonmode, phase)
    type (hits_type), INTENT(in) :: hits
    integer, INTENT(in) :: iFirstHit
    integer, INTENT(in) :: iLastHit
    logical, INTENT(in), dimension (kMaxHits,kMaxMult) :: isEdge
    type (fit_type), INTENT(in) :: fit
    type (uvfit_type), INTENT(in) :: uv
    logical, INTENT(in) :: muonmode
    type (zphifit_type), INTENT(out) :: phase
    ! Calls: NewZphiFit, addFit, tryPhase
end subroutine iteratePhase

private subroutine tryPhase (hits, hitsr, hitsp, hitsdp, r0, ifirst, ilast, ia1, im1, ia2, im2, iwind, minTant, maxTant, muonmode, phase)
    type (hits_type), INTENT(in) :: hits
    real (kind=R8), INTENT(in), dimension (:,:) :: hitsr
    real (kind=R8), INTENT(in), dimension (:,:) :: hitsp
    real (kind=R8), INTENT(in), dimension (:,:) :: hitsdp
    real (kind=R8), INTENT(in) :: r0
    integer, INTENT(in) :: ifirst
             first and last hits of subtrack
    integer, INTENT(in) :: ilast
             first and last hits of subtrack
    integer, INTENT(in) :: ia1
             first anchor
    integer, INTENT(in) :: im1
             first anchor
    integer, INTENT(in) :: ia2
             last anchor
    integer, INTENT(in) :: im2
             last anchor
    integer, INTENT(in) :: iwind
             anchor windings
    real (kind=R8), INTENT(in) :: minTant
             minimum and maximum tan(theta) allowed by cluster size
    real (kind=R8), INTENT(in) :: maxTant
             minimum and maximum tan(theta) allowed by cluster size
    logical, INTENT(in) :: muonmode
             true if particle is muon
    type (zphifit_type), INTENT(out) :: phase
    ! Calls: NewZPhiFit, fgkodie
end subroutine tryPhase

private function angleDiff (a, b)
    real (kind=R8), INTENT(in) :: a
    real (kind=R8), INTENT(in) :: b
    real (kind=R8) :: angleDiff
end function angleDiff

private function reduce (a)
    real (kind=R8), INTENT(in) :: a
    real (kind=R8) :: reduce
end function reduce

private function rad2deg (a)
    real (kind=R8), INTENT(in) :: a
    real (kind=R8) :: rad2deg
end function rad2deg

private function deg2rad (a)
    real (kind=R8), INTENT(in) :: a
    real (kind=R8) :: deg2rad
end function deg2rad

private function jsign (a) result (ival)
    real (kind=R8), INTENT(in) :: a
    integer :: ival
end function jsign

private function sqr (x) result (x2)
    real (kind=R8), INTENT(in) :: x
    real (kind=R8) :: x2
end function sqr

private subroutine Fgkodie ()
    ! Calls: kerror
end subroutine Fgkodie
 SUBROUTINE fgkodie

 "call fgkodie" prints the current run and event number and kills mofia.

 It is used like C "assert": to bomb out on internal errors
 and other unrecoverable or unhandled situations, where processing
 further events is dangerous and meaningless unless the asserted
 condition is fixed. K.O.

private function kcutrq2 (x) result (x2)
    real (kind=R8), INTENT(in) :: x
    real (kind=R8) :: x2
end function kcutrq2