Module Classify_mod

module Classify_mod

        ! Uses
    use Precision_mod
    use Chambers_mod
    use Det_Geom_mod
    use tdc_mod
    use unp_mod
    use unpmc_mod
    use Namelist_mod
    use Track_mod
    use Pattern_mod
    use Hists_mod
    use Skim_mod
    use trackrange_mod
    use Calibrations_mod, only: BField

        ! Variables
    integer (kind=I4), private, DIMENSION(2) :: nDCu
    integer (kind=I4), private, DIMENSION(2) :: nDCv
    integer (kind=I4), private, DIMENSION(2) :: nPairs
    integer (kind=i4), private, dimension (MaxWindows) :: MultiWindows
    integer (kind=i4), private :: iWinTracks
    integer (kind=i4), public, dimension (MaxWindows) :: nPlanesHit
    logical, public :: Selected
    real (kind=r4), public :: cp_time
    real (kind=r4), public :: rf_time
    real (kind=r4), public :: inavg
    real (kind=r4), public :: outavg
    real (kind=R4), public, dimension (MaxWindows) :: tdcwidavg
    real (kind=R4), public, dimension (MaxWindows) :: tdcwidavgwi
    real (kind=R4), public, dimension (MaxWindows) :: tdcwidavgpl
    logical, public, dimension (MaxWindows) :: TargetStop

        ! Subroutines and functions
    public function Classifying (Iteration) result (EventSelected)
    private subroutine EstimatePolarization ()
    private subroutine CalcTDCwidavg ()
    private subroutine FillDCTDCWidthHistos ()
    private subroutine EvalWindows ()
    private subroutine SumHitPlanes (iWindow)
    private subroutine EvalEvent ()
    public subroutine IdentifyPositrons ()
    public subroutine EvalTwoHits ()
    private subroutine EvalSpaceCharge ()
    public subroutine PostFitClassify (iStat)
    private subroutine FillFitFailHistMC (iFitFail)
    private subroutine FillMCEventTypeHist (iStat)
    public subroutine SelectTriggerParticle (UseTrigger)

end module Classify_mod
 MODULE Classify_mod
==============================================================================
 Author: Jim, Blair, Steven, Rob
 Initial Revision Feb. 27, 2002 ... Jim,Blair,Steven
 Last Updated: Aug. 2002 ... Rob
------------------------------------------------------------------------------
 This module contains subroutines used to classify events and identify 
 paritcles.  The main function of this module is to classify events before 
 running time expensive track fitting code (function classifying).  After 
 the event is tracked, there is also a post fit classification subroutine, to 
 decide what type of event has been fit (subroutine postFitClassify).  There 
 are also subroutines to evaluate the space charge effect (subroutine 
 EvalTwoHits), and to look do some studies of beam positrons (subroutine
 IdentifyPositrons).

Description of Variables

integer (kind=I4), private, DIMENSION(2) :: nDCu

integer (kind=I4), private, DIMENSION(2) :: nDCv

integer (kind=I4), private, DIMENSION(2) :: nPairs

integer (kind=i4), private, dimension (MaxWindows) :: MultiWindows
$  INTEGER(i4)::nDeltas

integer (kind=i4), private :: iWinTracks

integer (kind=i4), public, dimension (MaxWindows) :: nPlanesHit

logical, public :: Selected

real (kind=r4), public :: cp_time

real (kind=r4), public :: rf_time

real (kind=r4), public :: inavg

real (kind=r4), public :: outavg

real (kind=R4), public, dimension (MaxWindows) :: tdcwidavg

real (kind=R4), public, dimension (MaxWindows) :: tdcwidavgwi

real (kind=R4), public, dimension (MaxWindows) :: tdcwidavgpl

logical, public, dimension (MaxWindows) :: TargetStop

Description of Subroutines and Functions

public function Classifying (Iteration) result (EventSelected)
    integer (kind=I4), INTENT(in) :: Iteration
    logical :: EventSelected
    ! Calls: CalcTDCWidAvg, EstimatePolarization, EvalEvent, EvalWindows, FillDCTDCWidthhistos, FillMCEventTypeHist, FillSkim, HF1, hf1
end function Classifying
 FUNCTION Classifying(iteration) RESULT (EventSelected)
----------------------------------------------------------------------------
 This subroutine classifies events before tracking is done.  The classifying
 also classifies each window by calling EvalWindows, which fills the 
 window()%wtype values for each window.  The events are classified into the
 categories: Simple Clean (EventType==1), Time Clean (EventType==2), Simple
 Dirty (EventType==3), and ComplexDirty (EventType==4).  Right now, there
 are no simple dirty events, and simple clean + time clean + complex dirty =
 all events.  Also the Classifying code counts how many windows appear to 
 be beam positrons (sort of works), and how many deltas (not working yet)
 there are.  iteration = 1 for initial call of classifying.
 
 Classifying is recalled when delta removal or overlap handling has
 modified the windows.  iteration = 2 for second call.

 Once classifying is done, events can be selected for analysis, and or 
 skimming.  By default no skimming is done, and all events are analyzed.  To
 select certain event types to analyze, select how the classifying code 
 results will be used by setting name Classify SelectEvent to:
   SelectEvent Value   Outcome
   -----------------   -------
      -1               Analyze events which were NOT
                        selected by Classifying
       0 (default)     Classify events, but analyze
                        all events
       1               Analyze events which were
                        selected by Classifying
   any other value     Don't run classifying code, 
    (eg. 2)             analyze all events                      
 The type of event selected is set by setting one or more of the following 
 namelist values to be true:
   name Classify SelectSimpleClean  ! analyze only simple clean events     
   name Classify SelectTimeClean    ! analyze only time clean events       
   name Classify SelectOverlapDirty ! analyze only overlap events      
   name Classify SelectComplexDirty ! analyze only complex dirty events    
   name Classify SelectBeamPositron ! analyze events which contain beam e+ 
   name Classify SelectDelta        ! analyze events which contain deltas  

 Also, the type of events to skim can be independently selected by 
 setting all the skim namelists to set up skimming channels, and setting
 one of the following namelist variables to true:
   name skim SkimSimpleClean  ! skim only simple clean events to channel 1
   name skim SkimTimeClean    ! skim only time clean events to channel 2
   name skim SkimOverlapDirty ! skim only overlap events to channel 3
   name skim SkimComplexDirty ! skim only complex dirty events to channel 4
   name skim SkimBeamPositron ! skim events which contain beam e+ to channel 5
   name skim SkimDelta        ! skim events which contain deltas to channel 6

 This routine returns the a logical, which says whether the event was
 selected for analysis.

private subroutine EstimatePolarization ()
end subroutine EstimatePolarization

private subroutine CalcTDCwidavg ()
end subroutine CalcTDCwidavg
 Subroutine CalcTDCwidavg, by Blair Jamieson.
 In each window, sum all the hits, then divide by number
 of hits, number of wires, and number of planes, for three
 different averages.  These are useful for particle ID.
 -July 2002

private subroutine FillDCTDCWidthHistos ()
    ! Calls: hf1, hf2
end subroutine FillDCTDCWidthHistos
 Subroutine FillDCTDCWidthHistos, by Blair Jamieson.
 Fill histograms of the TDC width averages calculated in
 CalcTDCwidavg.  -July 2002

private subroutine EvalWindows ()
    ! Calls: HFF1, SumHitPlanes, hf1, hf2
end subroutine EvalWindows
 Subroutine EvalWindows
----------------------------------------------------------------------------
 EvalWindows classifies each window by particle type.
 Types (stored in WType(iWindow)) are:
            1     Muon
            2     UpStream Decay
            3     DownStream Decay
            4     Beam Positron
            5     Empty
            6     Overlap involved

            10    Trackable UpStream, a few DownStream Hits
            11    Trackable DownStream, a few UpStream Hits
            12    Trackable UpStream after "muon" and "decay"
            13    Trackable DownStream after "muon" and "decay"
            14    Trackable DownStream prior to muon
            15    Pass through the detector, but not beam positron
            16    DC clusters but no PC clusters.
 Type 12 could be the signal of a second muon. 
 Then 12 or 13 could signal a second decay.

 The value of TargetStop is also set for each window, regardless of
 whether the window is a muon or not.

private subroutine SumHitPlanes (iWindow)
    integer (kind=I4), INTENT(in) :: iWindow
end subroutine SumHitPlanes
----------------------------------------------------------------------------
 SumHitPlanes counts hits and hit pairs in the upstream and 
 downstream halves of the detector.  -Jim Musser.

 Rewritten version of SumHitPlanes which defines 'upstream' 
 'downstream' relative to triggerLastZ (from trackrange_mod),
 not the target.
 (Note that Classifying() sets triggerLastZ to 0 at the start
 if enableMuStopGuess is false.) -RPM 02.07.05

private subroutine EvalEvent ()
end subroutine EvalEvent
$  !!----------------------------------------------------------------------------
$  !! FindDelta steers the search for delta tracks and fills the structure Delta.
$  SUBROUTINE FindDelta(iWindow)
$
$    IMPLICIT NONE
$    INTEGER(I4), INTENT(in)::iWindow
$    INTEGER(I4)::iLine, iCl, iPair, jPair, nLines
$    INTEGER(I4), PARAMETER::MaxLines = 10, MinPairs = 5
$    LOGICAL::ClUsed
$    REAL(R8), PARAMETER::LineRes = 2.0
$!    TYPE::line_type
$!       INTEGER(I4)::FirstPair, LastPair, nPairs
$!       REAL(R8)::u, v
$!    END TYPE line_type
$    TYPE(delta_type), DIMENSION(MaxLines), TARGET::Line
$    TYPE(cluster_type), DIMENSION(:), POINTER::ClP
$
$    nLines = 0
$    jPair = 0
$    ClP => PCCl(iWindow, :)
$    CALL PairSearch(PPairLim(1, 1), PPairLim(1, 2))
$    ClP => DCCl(iWindow, :)
$    CALL PairSearch(DPairLim(1, 1), DPairLim(1, 3))
$    ClP => PCCl(iWindow, :)
$    CALL PairSearch(PPairLim(1, 2) + 1, PPairLim(2, 2) - 1)
$    ClP => DCCl(iWindow, :)
$    CALL PairSearch(DPairLim(2, 3), DPairLim(2, 1))
$    ClP => PCCl(iWindow, :)
$    CALL PairSearch(PPairLim(2, 2), PPairLim(2, 1))
$    ChLine: DO iLine = 1, nLines
$       IF(Line(iLine) % nPairs >= MinPairs .AND. &
$            REAL(Line(iLine) % nPairs)/&
$            REAL((Line(iLine) % LastPair - Line(iLine) % FirstPair + 1)) &
$            > 0.7) THEN
$          ! Delta
$          IF(nDeltas < MaxDeltas) THEN
$             nDeltas = nDeltas + 1
$             Delta(nDeltas) % iWindow = iWindow
$             Delta(nDeltas) % u = Line(iLine) % u
$             Delta(nDeltas) % v = Line(iLine) % v
$             IF(Line(iLine) % LastPair < 15) THEN ! UpStream delta
$                Delta(nDeltas) % Stream = 'U'
$             ELSE IF(Line(iLine) % FirstPair > 14) THEN ! DnStream delta
$                Delta(nDeltas) % Stream = 'D'
$             ELSE ! Both Streams delta
$                Delta(nDeltas) % Stream = 'B'
$             END IF
$          ELSE
$             EXIT ChLine
$          END IF
$       END IF
$    END DO ChLine
$
$  CONTAINS
$
$    !!--------------------------------------------------------------------------
$    !! PairSearch attempts to find clusters with a common (u, v) coordinate
$    !! indicative of a delta track.
$    SUBROUTINE PairSearch(FirstPair, LastPair)
$
$      IMPLICIT NONE
$      INTEGER(I4), INTENT(in)::FirstPair, LastPair
$      TYPE(coord_type), POINTER::CrdP
$      TYPE(delta_type), POINTER::LineP
$
$      ChPair: DO iPair = FirstPair, LastPair
$         jPair = jPair + 1
$         ChCluster: DO iCl = 1, ClP(iPair) % nCl
$            IF(ClP(iPair) % nCl == 1) CYCLE
$            CrdP => ClP(iPair) % Coord(iCl)
$            ClUsed = .FALSE.
$            ChTmpLine: DO iLine = 1, nLines
$               LineP => Line(iLine)
$               IF(ABS(CrdP % uBar - LineP % u) < LineRes &
$                    .AND. ABS(CrdP % vBar - LineP % v) < &
$                    LineRes) THEN
$                  ! Cluster is on line
$                  LineP % nPairs = LineP % nPairs + 1
$                  LineP % LastPair = jPair
$                  LineP % u = (LineP % u)*&
$                       (LineP % nPairs - 1)/LineP % nPairs + &
$                       (CrdP % uBar)/(LineP % nPairs)
$                  LineP % v = (LineP % v)*&
$                       (LineP % nPairs - 1)/LineP % nPairs + &
$                       (CrdP % vBar)/(LineP % nPairs)
$                  ClUsed = .TRUE.
$                  EXIT ChTmpLine
$               END IF
$            END DO ChTmpLine
$            IF(.NOT.ClUsed) THEN
$               ! Start new line
$               IF(nLines < MaxLines) THEN
$                  nLines = nLines + 1
$                  LineP => Line(nLines)
$                  LineP % u = CrdP % uBar
$                  LineP % v = CrdP % vBar
$                  LineP % FirstPair = jPair
$                  LineP % LastPair = jPair
$                  LineP % nPairs = 1
$               END IF
$            END IF
$         END DO ChCluster
$      END DO ChPair
$
$    END SUBROUTINE PairSearch
$    !---------------------------------------------------------------------------
$
$  END SUBROUTINE FindDelta

 SUBROUTINE EvalEvent
----------------------------------------------------------------------------
 EvalEvent classifies events according to the following:
 EventType =
             1(Simple Clean Event)    Muon and Decay, well separated in time
             2(Time Clean Event)      Muon, Decay, and Beam Positrons, well
                                           separated in time
             3(Overlap Dirty Event)   Two or more windows overlapped
             4(Complex Dirty Event)   Miscelaneous
             5(Time Clean, Non Muon Trigger)

public subroutine IdentifyPositrons ()
    ! Calls: HFF1, HFF2
end subroutine IdentifyPositrons
 SUBROUTINE IdentifyPositrons
------------------------------------------------------------------------------
 Identify Positrons looks up some variables which might help in identifying 
 positrons.  The variables of interest which are calculated by this subroutine
 are: (Blair 2002)
      cp_time  =  capacitive probe time for this event
      rf_time  =  rf time for this event
      nPlanesHit(iwindow) = number of planes hit by window
      iWinTracks = number of windows with > 30 planes hit
      MultiWindows(iWinTracks) = window numbers with > 30 planes hit
      inavg    = average pc time for incident straight through particle
      outavg   = average pc time for outgoing straight through particle

public subroutine EvalTwoHits ()
    ! Calls: EvalSpaceCharge, HFF1, HFF2
end subroutine EvalTwoHits
 SUBROUTINE SkimMuonHits
------------------------------------------------------------------------------
 SkimMuonHits produces histograms designed to answer question concerning
 whether or not wires hit by the muon are dead for an extended period of time.
 This subroutine was written and executed on the evening of Dec. 17, 2001 for
 the benefit of the NSERC review.  The histograms showed no noticeable
 suppression of hits due to overlap of muon and positron hits.  In fact,
 using the standard cut on prompt decays at 1 us, the code showed the
 lifetime of the muon to be 2.19 us.
 Authors: Blair Jamieson, Rob MacDonald, Jim Musser

 EvalTwoHits is based on the subroutine SkimMuonHits which is described above.
 Additions to EvalTwoHits include using firstguess and tracking information to
 estimate the separation of the two hits on a wire. (Blair 2002)

private subroutine EvalSpaceCharge ()
    ! Calls: HF1, HF2, hf1, hf2
end subroutine EvalSpaceCharge
 SUBROUTINE EvalSpaceCharge
------------------------------------------------------------------------------
 EvalSpaceCharge produces histograms designed to answer question concerning
 whether or not wires hit by the muon are dead for an extended period of time.

 The algorithm is:
   loop over windows looking for window categorized as a beam positron
      if the beam positron was successfully fit then
         loop through hits in muon window
            if muon hit is in the same cell as the positron passes through
               find the distance from the muon hit to the positron path
               find path length of positron through cell
               look for a positron hit in the same cell as the muon 
               find the distance between the two hits if there are two hits
               fill histograms

public subroutine PostFitClassify (iStat)
    integer (kind=I4), INTENT(in) :: iStat
    ! Calls: FillFitFailHistMC, HFF1
end subroutine PostFitClassify
 SUBROUTINE PostFitClassify
------------------------------------------------------------------------------
 Right now this isn't too useful. It tries to classify fits based on how
 good the first guess and helix fits are.
 

private subroutine FillFitFailHistMC (iFitFail)
    integer (kind=I4), INTENT(in) :: iFitFail
    ! Calls: HFF1, HFF2
end subroutine FillFitFailHistMC

private subroutine FillMCEventTypeHist (iStat)
    integer (kind=I4), INTENT(in) :: iStat
    ! Calls: HF1, HF2
end subroutine FillMCEventTypeHist

public subroutine SelectTriggerParticle (UseTrigger)
    logical, INTENT(OUT) :: UseTrigger
    ! Calls: HF2
end subroutine SelectTriggerParticle
 SUBROUTINE SelectTriggerParticle( UseTrigger )
------------------------------------------------------
 This subroutine decides based on CPTOF, and MU SCINTILLATOR 
 M12 width whether to process rest of this event.  The 
 return value is TRUE is the event is inside the selected
 gates in TCAP, and M12 width, and false otherwise.

 namelists to set are:
  name classify selectTCAPmin = ...
  name classify selectTCAPmax = ...
  name classify selectM12min = ...
  name classify selectM12max = ...