*--------1---------2---------3---------4---------5---------6---------7--

* Distribution of the SPACAL clusters with
* respect to the BST tracks to study event
* signatures and background patterns

* declare variables
*----------------------------------
      real function DIFF(Energy)

      include './Include/common.inc'
      include './Include/padgeom.inc'

      real Energy
      external CTVDCA, CNXY
      real RIN(3), VTX(3), ROUT(5)

      real XH, YH, ZH, RH, PH, TH, RS, PS, PT, Pi
      real RO(4), PO(4), RE(4), PE(4)
      real LLIM, ULIM, HPAR

      integer NHITE(6), NHITO(6), Sec_Trk(12), NSEC
      integer NUM(12, 8), DEN(12, 8)

      logical CAND, TRACK

* initialize variables
*----------------------------------
      parameter (Pi = 3.14159265358979)

      CAND = .FALSE.
      TRACK = .FALSE.

      NSEC = 0
      EFFECT = 1

      do I  = 1, 6
        NHITE(I) = 0
	NHITO(I) = 0
      end do

      do I = 1, 12
        Sec_Trk(I) = 0
      end do


* book histograms
*----------------------------------
      if (IDNEVT.eq.1) then
        call HBOOK1(1, 'F difference', 20,0.,360.,0)
        call HBOOK1(2, 'R difference', 20,0.,35.,0)

        call HBOOK1(3, 'F difference', 20,0.,360.,0)
        call HBOOK1(4, 'R difference', 20,0.,35.,0)

        call HBOOK1(5, 'F difference', 20,0.,360.,0)
        call HBOOK1(6, 'R difference', 20,0.,35.,0)

        call HBOOK1(7, 'F difference', 20,0.,360.,0)
        call HBOOK1(8, 'R difference', 20,0.,35.,0)
      end if

* strip track in the pad acceptance
*----------------------------------
      RIN(1) = INRAD(1)
      RIN(2) = PHBST(1)
      RIN(3) = THBST(1)
      
      VTX(1) = XVBST(1)
      VTX(2) = YVBST(1)
      VTX(3) = ZVBST(1)

      call CTVDCA(RIN, VTX, ROUT)

      do J = 1, 4

* validating even sector
*----------------------------------
        ZH = (-1.) * EvenZPos(J) / 10.

        call CNXY (ROUT, ZH, XH, YH)

        RH = 10. * sqrt(XH**2 + YH**2)
        PH = 180. * atan2(YH, XH) / Pi
        if (PH.lt.0.) PH = PH + 360.

        do K = 1, 6
          do L = 1, 8
            if ((RH.gt.RadDiv(L)).and.(RH.lt.RadDiv(L + 1)).and.
     +          (PH.gt.EvenMinPhi(K)).and.(PH.lt.EvenMaxPhi(K)))
     +      then
              NHITE(K) = NHITE(K) + 1
              RE(J) = RH
              PE(J) = PH
            end if
          end do
        end do


        ZH = (-1.) * OddZPos(J) / 10.

        call CNXY (ROUT, ZH, XH, YH)

        RH = 10. * sqrt(XH**2 + YH**2)
        PH = 180. * atan2(YH, XH) / Pi
        if (PH.lt.0.) PH = PH + 360.

* validating odd sector
*----------------------------------
        do K = 1, 6
          do L = 1, 8
            if ((RH.gt.RadDiv(L)).and.(RH.lt.RadDiv(L + 1)).and.
     +          (PH.gt.OddMinPhi(K)).and.(PH.lt.OddMaxPhi(K)))
     +      then
              NHITO(K) = NHITO(K) + 1
              RO(J) = RH
              PO(J) = PH
            end if
          end do
        end do

      end do

* track candidate
*----------------------------------
      do I = 1, 6
        J = 2 * I - 1
	K = 2 * I
	
	if (NHITE(I).eq.4) Sec_Trk(J) = 1
	if (NHITO(I).eq.4) Sec_Trk(K) = 1
      end do

      do I = 1, 12
        if (Sec_Trk(I).eq.1) NSEC = NSEC + 1
      end do

      if (NSEC.eq.1) CAND = .TRUE.

* check the SPACAL energy
*----------------------------------      
      TH = atan2(RS, ZH)
      PT = 0.3 / 100. / INRAD(1) / sin(TH)

      if (NSEL.le.0.or.ESEL(1).lt.Energy) CAND = .FALSE.

* link BST tracks to SPACAL clusters
*----------------------------------
      ZH = (-1.) * SpacalPos / 10.

      call CNXY (ROUT, ZH, XH, YH)

      RH = sqrt(XH**2 + YH**2)
      PH = 180. * atan2(YH, XH) / Pi
      if (PH.lt.0.) PH = PH + 360.

      RS = sqrt(XSEL(1)**2 + YSEL(1)**2)
      PS = 180. * atan2(YSEL(1), XSEL(1)) / Pi
      if (PS.lt.0.) PS = PS + 360.


* filling histograms
*----------------------------------
      HPAR = abs(PH-PS)
      call HF1(1, HPAR, 1.)
      if (abs(PH-PS).le.22.) call HF1(3, HPAR, 1.)
      if (abs(RH-RS).le.4.)  call HF1(5, HPAR, 1.)
      if ((abs(PH-PS).le.22.).and.(abs(RH-RS).le.4.))
     +     call HF1(7, HPAR, 1.)

      HPAR = abs(RH-RS)
      call HF1(2, HPAR, 1.)
      if (abs(RH-RS).le.4.) call HF1(4, HPAR, 1.)
      if (abs(PH-PS).le.22.) call HF1(6, HPAR, 1.)
      if ((abs(PH-PS).le.22.).and.(abs(RH-RS).le.4.))
     +     call HF1(8, HPAR, 1.)

      return
      end