*file member=acf_WindowCut  library=snoman language=fortran77 date=05:Jul:2001
 <<<      subroutine acf_WindowCut
         +         (vtxPos,nHitPMT,PMTtime,PMTlist,FirstRTCPMT,LastRTCPMT)
     
    *     acf_WindowCut: Sorts the PMTs in ascending order by residual hit time,
    *       and finds the largest subset of the PMTs that have hit times within
    *       the PMT time window of each other.
     
    *     Contact:  W. Heintzelman, Penn
    *     Verified: 
    *     Refereed: 
     
    *     Parameters:-
    *     ==========
    *     vtxPos    in      Position and time of event
    *     nHitPMT   in      Number of hit PMTs
    *     PMTtime   in/out  Times of PMT hits (returned sorted by residual time)
    *     PMTlist   in/out  List of hit PMTs (returned sorted by residual time)
    *     FirstRTCPMT out   Index of first PMT in the residual time cut window
    *     LastRTCPMT  out   Index of last  PMT in the residual time cut window
    *     DISCARD   in      = .FALSE.
    *               out     Set = .TRUE. if PMTpos array could not be filled
    *         
     
    *     Common Block Access:-
    *     ===================
    *        , Banks: 
     
    *     Specification:-
    *     =============
    *       Sort the PMTs in ascending order by residual hit time.
    *       Find the largest subset of the PMTs that have residual hit times
    *       within the PMT time window of each other, and save the indices of
    *       the first and last PMTs in the subset.
     
    *     Program Notes:-
    *     =============
     
    *     Revision History:-
    *     ================
    *     4.02   W. Heintzelman   First version
    *            W. Heintzelman   Second version
    *     5.00   J. Formaggio     Expanded pmt numbers to accomodate for muon hits (>1000)
     
          IMPLICIT NONE
     
          include 'zunit.inc'
          include 'acf_evalconst.inc'
          include 'ge_mnemonics.inc'
          include 'ge_pmts.inc'
     
    *     Argument Declarations:-
    *     =====================
            integer nHitPMT,PMTlist(nHitPMT),FirstRTCPMT,LastRTCPMT
            real vtxPos(4),PMTtime(nHitPMT)
     
    *ENDHEADER
     
    *     Local Variable Declarations:-
    *     ===========================
            integer maxPMT
            parameter (maxPMT=10000)
     
            integer itag(maxPMT),i,j,ninclmax,nincl,itemp(maxPMT)
            real temp(maxPMT),rtemp(maxPMT),Tresid(maxPMT),dsq
     
            if(nHitPMT.gt.maxPMT) then
                write(iqlog,*)
         +         'Fatal error in acf_WindowCut.  Input number of PMT hits'
         +              ,' exceeds dimension provided.  nHitPMT =',nHitPMT
                stop
            endif
     
            ! Calculate the PMT hit time residuals
            ! (21.75 is the velocity of light in water recommended by Josh in
            !  cm/ns.
            do i=1,nHitPMT
                dsq = 0.
                do j=1,3
                    ! dsq = dsq +(PMTcoord(j,PMTlist(i))-vtxPos(j))**2
                    dsq = dsq +(RGE_PMT_XYZ(j,PMTlist(i))-vtxPos(j))**2
                enddo
                Tresid(i) = PMTtime(i)-vtxPos(4)-sqrt(dsq)/21.75
            enddo
     
    *        write(iqlog,*) 'PMT list at entry:'
    *        do i=1,nHitPMT
    *            write(iqlog,'(i3,2f10.2,i6)')
    *     +               i,Tresid(i),PMTtime(i),PMTlist(i)
    *        enddo
     
            ! call fsort(Tresid,itag,nHitPMT)
            ! Call to my routine replaced by call to CERN library routine:
            call sortzv(Tresid,itag,nHitPMT,1,0,0)    ! Sort the PMT hits by
                                                      ! residual hit time
            ! Save PMT times and IDs in temporary storage:
            do i=1,nHitPMT
                temp(i) = PMTtime(i)
                rtemp(i) = Tresid(i)
                itemp(i) = PMTlist(i)
            enddo
            ! Now put the PMTs in ascending order by residual hit time:
            do i=1,nHitPMT
                PMTtime(i) = temp(itag(i))
                Tresid(i)  = rtemp(itag(i))
                PMTlist(i) = itemp(itag(i))
            enddo
     
    *        write(iqlog,*)
    *        write(iqlog,*) 'acf_Windowcut, PMT list after sort:'
    *        do i=1,nHitPMT
    *             write(iqlog,'(i3,2f10.2,i6)')
    *     +                i,Tresid(i),PMTtime(i),PMTlist(i)
    *         enddo
     
    *       Spread of 8 ns for the residuals defined in acf_evalconst.inc as
    *       parameter TWindow.
            ninclmax=0
            FirstRTCPMT=0
            do i=1,nHitPMT  ! Loop over possible earliest hits in the subset
                nincl=0
                ! Loop over later hits, counting the number within TWindow of the
                ! ith hit:
                do j=i,nHitPMT
                    if( (Tresid(j)-Tresid(i)) .le. TWindow ) then
                        nincl = nincl + 1
                    else
                        ! j=nHitPMT           ! Jump out of loop
                        goto 100
                    endif
                enddo
      100       continue
     
                ! Save the index to the ith hit if there are more hits within
                ! TWindow of it than for previous cases:
                if(nincl.gt.ninclmax) then
                    ninclmax = nincl        ! Number of PMTs in the window
                    FirstRTCPMT=i           ! First one included
                endif
            enddo
     
            LastRTCPMT = FirstRTCPMT -1 + ninclmax
     
            ! Debug printout
    *        write(iqlog,*) 'acf_WindowCut: FirstRTCPMT,LastRTCPMT =',
    *     +          FirstRTCPMT,LastRTCPMT
    *        do i=1,nHitPMT
    *          write(iqlog,*) i,tresid(i),pmttime(i),
    *     +                   (RGE_PMT_XYZ(j,PMTlist(i)),j=1,3)
    *        enddo
     
            return
            end
    *endfile member=acf_WindowCut