*****************************************************************************
      INTEGER FUNCTION TOFTRK()
c     
      IMPLICIT NONE
C     
      include  'input_tof.txt'
      include  'output_tof.txt'
      include  'tofcomm.txt'
C     
      
c     =======================================
c     variables for tracking routine
c     =======================================
      integer NPOINT_MAX
      parameter(NPOINT_MAX=100)

c     define TOF Z-coordinates
      integer NPTOF
      parameter (NPTOF=6)
      DOUBLE PRECISION ZTOF(NPTOF)
      DATA ZTOF/53.74,53.04,23.94,23.44,-23.49,-24.34/ !Sergio 9.05.2006

      integer itof

      DOUBLE PRECISION al_p(5),
     &     xout(NPOINT_MAX),yout(NPOINT_MAX),zin(NPTOF)

      INTEGER IFAIL
      REAL dx,dy,dr,ds
      REAL yhelp,xdummy,xkorr0,xhelp,xhelp1,xhelp2
      REAL c1,c2,sw,sxw,w_i
      INTEGER icount

      REAL beta_mean


      REAL theta12,theta13,theta23
C--   DATA ZTOF/53.74,53.04,23.94,23.44,-23.49,-24.34/ !Sergio 9.05.2006
      REAL tofarm12
      PARAMETER (tofarm12 = 29.70) ! from 53.39 to 23.69
      REAL tofarm23	   
      PARAMETER (tofarm23 = 47.61)  ! from 23.69 to -23.92
      REAL tofarm13
      PARAMETER (tofarm13 = 77.31)  ! from 53.39 to -23.92

      
      
      
      INTEGER ihelp
      REAL xkorr

      REAL yl,yh,xl,xh
C     
      REAL hmemor(9000000)
      INTEGER Iquest(100)
C     
      COMMON / pawcd / hmemor
      save / pawcd /
C     
      Common / QUESTd / Iquest
      save / questd /
C     
C     Begin !
C     
      TOFTRK = 0

*******************************************************************

      offset = 1
      slope = 2
      left = 1
      right = 2
      none_ev = 0
      none_find = 0
      tdc_ev = 1
      adc_ev = 1
      itdc = 1
      iadc = 2

      do i=1,13
         beta_a(i) = 100.
      enddo

C------ read  tracking routine
*     igoodevent = igoodevent+1
*     assigned input  parameters for track routine
*     1) Z-coordinates where the trajectory is evaluated
      do itof=1,NPTOF
         ZIN(itof) = ZTOF(itof)
      enddo
*     2) track status vector
C     COPY THE ALFA VECTOR FROM AL_PP TO AL_P FOR THE TRACK "T"
      do i=1,5
         AL_P(i) = al_pp(i)
      enddo
      if (al_p(5).eq.0.) THEN
         PRINT *,' TOF - WARNING F77: track with R = 0, discarded'
         GOTO 969
      ENDIF
*     -------- *** tracking routine *** --------
      IFAIL = 0
      call DOTRACK(NPTOF,ZIN,XOUT,YOUT,AL_P,IFAIL)
      if(IFAIL.ne.0)then
         print *,' TOF - WARNING F77: tracking failed '
         goto 969
      endif
*     ------------------------------------------

969   continue

C-------------------------------------------------------------
C-------  check which paddle penetrated the track  -----------
C-------------------------------------------------------------
c middle y (or x) position of the upper and middle ToF-Paddle
c       DATA tof11_x/ -17.85,-12.75,-7.65,-2.55,2.55,7.65,12.75,17.85/
c       DATA tof12_y/ -13.75,-8.25,-2.75,2.75,8.25,13.75/
c       DATA tof21_y/ -3.75,3.75/
c       DATA tof22_x/ -4.5,4.5/
c       DATA tof31_x/ -6.0,0.,6.0/
c       DATA tof32_y/ -5.0,0.0,5.0/
c
c  S11 8 paddles  33.0 x 5.1 cm
c  S12 6 paddles  40.8 x 5.5 cm
c  S21 2 paddles  18.0 x 7.5 cm
c  S22 2 paddles  15.0 x 9.0 cm
c  S31 3 paddles  15.0 x 6.0 cm
c  S32 3 paddles  18.0 x 5.0 cm


C--------------     S11 --------------------------------------

      tof11_i = none_find

      yl = -33.0/2.
      yh =  33.0/2.

      if ((yout(1).gt.yl).and.(yout(1).lt.yh)) then
      do i=1,8
      xl = tof11_x(i) - 5.1/2.
      xh = tof11_x(i) + 5.1/2.
      if ((xout(1).gt.xl).and.(xout(1).le.xh)) then
      tof11_i=i
      endif
      enddo
      endif

C--------------     S12 --------------------------------------

      tof12_i = none_find
 
      xl = -40.8/2.
      xh =  40.8/2.

      if ((xout(2).gt.xl).and.(xout(2).lt.xh)) then
      do i=1,6
      yl = tof12_y(i) - 5.5/2.
      yh = tof12_y(i) + 5.5/2.
      if ((yout(2).gt.yl).and.(yout(2).le.yh)) then
      tof12_i=i
      endif
      enddo
      endif

C--------------     S21 --------------------------------------

      tof21_i = none_find
 
      xl = -18./2.
      xh =  18./2.

      if ((xout(3).gt.xl).and.(xout(3).lt.xh)) then
      do i=1,2
      yl = tof21_y(i) - 7.5/2.
      yh = tof21_y(i) + 7.5/2.
      if ((yout(3).gt.yl).and.(yout(3).le.yh)) then
      tof21_i=i
      endif
      enddo
      endif

C--------------     S22 --------------------------------------

      tof22_i = none_find
 
      yl = -15./2.
      yh =  15./2.

      if ((yout(4).gt.yl).and.(yout(4).lt.yh)) then
      do i=1,2
      xl = tof22_x(i) - 9.0/2.
      xh = tof22_x(i) + 9.0/2.
      if ((xout(4).gt.xl).and.(xout(4).le.xh)) then
      tof22_i=i
      endif
      enddo
      endif

C--------------     S31 --------------------------------------

      tof31_i = none_find
 
      yl = -15.0/2.
      yh =  15.0/2.

      if ((yout(5).gt.yl).and.(yout(5).lt.yh)) then
      do i=1,3
      xl = tof31_x(i) - 6.0/2.
      xh = tof31_x(i) + 6.0/2.
      if ((xout(5).gt.xl).and.(xout(5).le.xh)) then
      tof31_i=i
      endif
      enddo
      endif

C--------------     S32 --------------------------------------

      tof32_i = none_find
 
      xl = -18.0/2.
      xh =  18.0/2.

      if ((xout(6).gt.xl).and.(xout(6).lt.xh)) then
      do i=1,3
      yl = tof32_y(i) - 5.0/2.
      yh = tof32_y(i) + 5.0/2.
      if ((yout(6).gt.yl).and.(yout(6).le.yh)) then
      tof32_i=i
      endif
      enddo
      endif

C----------------------------------------------------------------------
C--------------------Corrections on ADC-data -------------------------
C---------------------zenith angle theta  ---------------------------
C----------------------------------------------------------------------


      dx=0.
      dy=0.
      dr=0.
      theta13 = 0.

      if (xout(1).lt.100.) then
         dx = xout(1)-xout(6)
         dy = yout(1)-yout(6)
         dr = sqrt(dx*dx+dy*dy)
         theta13 = atan(dr/tofarm13)
      endif


C----------------------------------------------------------------------
C------------------angle and ADC(x) correction
C----------------------------------------------------------------------
C-----------------------------S1 --------------------------------

      yhelp=yout(1)

      IF (tof11_i.GT.none_find.AND.yhelp.lt.100) THEN

         i = tof11_i
         xdummy=tof11(left,i,iadc) 
         tof11(left,i,iadc) = tof11(left,i,iadc)*cos(theta13)
         if (tof11(left,i,iadc).lt.1000) then
            xkorr=adcx11(left,i,1)*exp(-yhelp/adcx11(left,i,2))
            xkorr0=adcx11(left,i,1)
            adc_c(ch11a(i),hb11a(i))=tof11(left,i,iadc)/xkorr
         endif

         tof11(right,i,iadc) = tof11(right,i,iadc)*cos(theta13)
         if (tof11(right,i,iadc).lt.1000) then
            xkorr=adcx11(right,i,1)*exp(yhelp/adcx11(right,i,2))
            xkorr0=adcx11(right,i,1)
            adc_c(ch11b(i),hb11b(i))=tof11(right,i,iadc)/xkorr
         endif
      ENDIF


      xhelp=xout(2)
      IF (tof12_i.GT.none_find.AND.xhelp.lt.100) THEN

         i = tof12_i
         tof12(left,i,iadc) = tof12(left,i,iadc)*cos(theta13)
         if (tof12(left,i,iadc).lt.1000) then
            xkorr=adcx12(left,i,1)*exp(-xhelp/adcx12(left,i,2))
            xkorr0=adcx12(left,i,1)
            adc_c(ch12a(i),hb12a(i))=tof12(left,i,iadc)/xkorr
         endif

         tof12(right,i,iadc) = tof12(right,i,iadc)*cos(theta13)
         if (tof12(right,i,iadc).lt.1000) then
            xkorr=adcx12(right,i,1)*exp(xhelp/adcx12(right,i,2))
            xkorr0=adcx12(right,i,1)
            adc_c(ch12b(i),hb12b(i))=tof12(right,i,iadc)/xkorr
         endif
      ENDIF

C-----------------------------S2 --------------------------------

      xhelp=xout(3)
      IF (tof21_i.GT.none_find.AND.xhelp.lt.100) THEN

         i = tof21_i
         tof21(left,i,iadc) = tof21(left,i,iadc)*cos(theta13)
         if (tof21(left,i,iadc).lt.1000) then
            xkorr=adcx21(left,i,1)*exp(-xhelp/adcx21(left,i,2))
            xkorr0=adcx21(left,i,1)
            adc_c(ch21a(i),hb21a(i))=tof21(left,i,iadc)/xkorr
         endif

         tof21(right,i,iadc) = tof21(right,i,iadc)*cos(theta13)
         if (tof21(right,i,iadc).lt.1000) then
            xkorr=adcx21(right,i,1)*exp(xhelp/adcx21(right,i,2))
            xkorr0=adcx21(right,i,1)
            adc_c(ch21b(i),hb21b(i))=tof21(right,i,iadc)/xkorr
         endif
      ENDIF

      yhelp=yout(4)
      IF (tof22_i.GT.none_find.AND.yhelp.lt.100) THEN

         i = tof22_i
         tof22(left,i,iadc) = tof22(left,i,iadc)*cos(theta13)
         if (tof22(left,i,iadc).lt.1000) then
            xkorr=adcx22(left,i,1)*exp(-yhelp/adcx22(left,i,2))
            xkorr0=adcx22(left,i,1)
            adc_c(ch22a(i),hb22a(i))=tof22(left,i,iadc)/xkorr
         endif

         tof22(right,i,iadc) = tof22(right,i,iadc)*cos(theta13)
         if (tof22(right,i,iadc).lt.1000) then
            xkorr=adcx22(right,i,1)*exp(yhelp/adcx22(right,i,2))
            xkorr0=adcx22(right,i,1)
            adc_c(ch22b(i),hb22b(i))=tof22(right,i,iadc)/xkorr
         endif
      ENDIF

C-----------------------------S3 --------------------------------

      yhelp=yout(5)
      IF (tof31_i.GT.none_find.AND.yhelp.lt.100) THEN

         i = tof31_i
         tof31(left,i,iadc) = tof31(left,i,iadc)*cos(theta13)
         if (tof31(left,i,iadc).lt.1000) then
            xkorr=adcx31(left,i,1)*exp(-yhelp/adcx31(left,i,2))
            xkorr0=adcx31(left,i,1)
            adc_c(ch31a(i),hb31a(i))=tof31(left,i,iadc)/xkorr
         endif

         tof31(right,i,iadc) = tof31(right,i,iadc)*cos(theta13)
         if (tof31(right,i,iadc).lt.1000) then
            xkorr=adcx31(right,i,1)*exp(yhelp/adcx31(right,i,2))
            xkorr0=adcx31(right,i,1)
            adc_c(ch31b(i),hb31b(i))=tof31(right,i,iadc)/xkorr
         endif
      ENDIF

      xhelp=xout(6)
      IF (tof32_i.GT.none_find.AND.xhelp.lt.100) THEN

         i = tof32_i
         tof32(left,i,iadc) = tof32(left,i,iadc)*cos(theta13)
         if (tof32(left,i,iadc).lt.1000) then
            xkorr=adcx32(left,i,1)*exp(-xhelp/adcx32(left,i,2))
            xkorr0=adcx32(left,i,1)
            adc_c(ch32a(i),hb32a(i))=tof32(left,i,iadc)/xkorr
         endif

         tof32(right,i,iadc) = tof32(right,i,iadc)*cos(theta13)
         if (tof32(right,i,iadc).lt.1000) then
            xkorr=adcx32(right,i,1)*exp(xhelp/adcx32(right,i,2))
            xkorr0=adcx32(right,i,1)
            adc_c(ch32b(i),hb32b(i))=tof32(right,i,iadc)/xkorr
         endif
      ENDIF

C-----------------------------------------------------------------------
C----------------------calculate Beta  ------------------------------
C-----------------------------------------------------------------------
C-------------------difference of sums  ---------------------------
C
C     DS = (t1+t2) - t3+t4)
C     DS = c1 + c2/beta*cos(theta)
C     c2 = 2d/c   gives c2 = 2d/(c*TDCresolution)  TDC=50ps/channel
C     =>  c2 =  ca.60 for  0.45 m    c2 = ca.109  for 0.81 m
C     since TDC resolution varies slightly c2 has to be calibrated

C     S11 - S31
      IF (tof11_i.GT.none_find.AND.tof31_i.GT.none_find) THEN
         xhelp1 = tof11(1,tof11_i,itdc)+tof11(2,tof11_i,itdc)
         xhelp2 = tof31(1,tof31_i,itdc)+tof31(2,tof31_i,itdc)
         ds = xhelp1-xhelp2
         ihelp=(tof11_i-1)*3+tof31_i
         c1 = k_S11S31(1,ihelp)
         c2 = k_S11S31(2,ihelp)         
	 if ((xhelp1.lt.8000.).and.(xhelp2.lt.8000)) 
     &   beta_a(1) = c2/(cos(theta13)*(ds-c1))
      ENDIF
      
C     S11 - S32
      IF (tof11_i.GT.none_find.AND.tof32_i.GT.none_find) THEN
         xhelp1 = tof11(1,tof11_i,itdc)+tof11(2,tof11_i,itdc)
         xhelp2 = tof32(1,tof32_i,itdc)+tof32(2,tof32_i,itdc)
         ds = xhelp1-xhelp2
         ihelp=(tof11_i-1)*3+tof32_i
         c1 = k_S11S32(1,ihelp)
         c2 = k_S11S32(2,ihelp)         
	 if ((xhelp1.lt.8000.).and.(xhelp2.lt.8000)) 
     &   beta_a(2) = c2/(cos(theta13)*(ds-c1))
      ENDIF

C     S12 - S31
      IF (tof12_i.GT.none_find.AND.tof31_i.GT.none_find) THEN
         xhelp1 = tof12(1,tof12_i,itdc)+tof12(2,tof12_i,itdc)
         xhelp2 = tof31(1,tof31_i,itdc)+tof31(2,tof31_i,itdc)
         ds = xhelp1-xhelp2
         ihelp=(tof12_i-1)*3+tof31_i
         c1 = k_S12S31(1,ihelp)
         c2 = k_S12S31(2,ihelp)
	 if ((xhelp1.lt.8000.).and.(xhelp2.lt.8000)) 
     &   beta_a(3) = c2/(cos(theta13)*(ds-c1))
      ENDIF

C     S12 - S32
      IF (tof12_i.GT.none_find.AND.tof32_i.GT.none_find) THEN
         xhelp1 = tof12(1,tof12_i,itdc)+tof12(2,tof12_i,itdc)
         xhelp2 = tof32(1,tof32_i,itdc)+tof32(2,tof32_i,itdc)
         ds = xhelp1-xhelp2
         ihelp=(tof12_i-1)*3+tof32_i
         c1 = k_S12S32(1,ihelp)
         c2 = k_S12S32(2,ihelp)
	 if ((xhelp1.lt.8000.).and.(xhelp2.lt.8000)) 
     &   beta_a(4) = c2/(cos(theta13)*(ds-c1))
      ENDIF

C     S21 - S31
      IF (tof21_i.GT.none_find.AND.tof31_i.GT.none_find) THEN
         xhelp1 = tof21(1,tof21_i,itdc)+tof21(2,tof21_i,itdc)
         xhelp2 = tof31(1,tof31_i,itdc)+tof31(2,tof31_i,itdc)
         ds = xhelp1-xhelp2
         ihelp=(tof21_i-1)*3+tof31_i
         c1 = k_S21S31(1,ihelp)
         c2 = k_S21S31(2,ihelp)
	 if ((xhelp1.lt.8000.).and.(xhelp2.lt.8000)) 
     &   beta_a(5) = c2/(cos(theta23)*(ds-c1))
      ENDIF

C     S21 - S32
      IF (tof21_i.GT.none_find.AND.tof32_i.GT.none_find) THEN
         xhelp1 = tof21(1,tof21_i,itdc)+tof21(2,tof21_i,itdc)
         xhelp2 = tof32(1,tof32_i,itdc)+tof32(2,tof32_i,itdc)
         ds = xhelp1-xhelp2
         ihelp=(tof21_i-1)*3+tof32_i
         c1 = k_S21S32(1,ihelp)
         c2 = k_S21S32(2,ihelp)
	 if ((xhelp1.lt.8000.).and.(xhelp2.lt.8000)) 
     &   beta_a(6) = c2/(cos(theta23)*(ds-c1))
      ENDIF

C     S22 - S31
      IF (tof22_i.GT.none_find.AND.tof31_i.GT.none_find) THEN
         xhelp1 = tof22(1,tof22_i,itdc)+tof22(2,tof22_i,itdc)
         xhelp2 = tof31(1,tof31_i,itdc)+tof31(2,tof31_i,itdc)
         ds = xhelp1-xhelp2
         ihelp=(tof22_i-1)*3+tof31_i
         c1 = k_S22S31(1,ihelp)
         c2 = k_S22S31(2,ihelp)
	 if ((xhelp1.lt.8000.).and.(xhelp2.lt.8000)) 
     &   beta_a(7) = c2/(cos(theta13)*(ds-c1))
      ENDIF
      
C     S22 - S32
      IF (tof22_i.GT.none_find.AND.tof32_i.GT.none_find) THEN
         xhelp1 = tof22(1,tof22_i,itdc)+tof22(2,tof22_i,itdc)
         xhelp2 = tof32(1,tof32_i,itdc)+tof32(2,tof32_i,itdc)
         ds = xhelp1-xhelp2
         ihelp=(tof22_i-1)*3+tof32_i
         c1 = k_S22S32(1,ihelp)
         c2 = k_S22S32(2,ihelp)
	 if ((xhelp1.lt.8000.).and.(xhelp2.lt.8000)) 
     &   beta_a(8) = c2/(cos(theta13)*(ds-c1))
      ENDIF

C     S11 - S21
      IF (tof11_i.GT.none_find.AND.tof21_i.GT.none_find) THEN
         xhelp1 = tof11(1,tof11_i,itdc)+tof11(2,tof11_i,itdc)
         xhelp2 = tof21(1,tof21_i,itdc)+tof21(2,tof21_i,itdc)
         ds = xhelp1-xhelp2
         ihelp=(tof11_i-1)*2+tof21_i
         c1 = k_S11S21(1,ihelp)
         c2 = k_S11S21(2,ihelp)
	 if ((xhelp1.lt.8000.).and.(xhelp2.lt.8000)) 
     &   beta_a(9) = c2/(cos(theta12)*(ds-c1))
      ENDIF
      
C     S11 - S22
      IF (tof11_i.GT.none_find.AND.tof22_i.GT.none_find) THEN
         xhelp1 = tof11(1,tof11_i,itdc)+tof11(2,tof11_i,itdc)
         xhelp2 = tof22(1,tof22_i,itdc)+tof22(2,tof22_i,itdc)
         ds = xhelp1-xhelp2
         ihelp=(tof11_i-1)*2+tof22_i
         c1 = k_S11S22(1,ihelp)
         c2 = k_S11S22(2,ihelp)
	 if ((xhelp1.lt.8000.).and.(xhelp2.lt.8000)) 
     &   beta_a(10) = c2/(cos(theta12)*(ds-c1))
      ENDIF

C     S12 - S21
      IF (tof12_i.GT.none_find.AND.tof21_i.GT.none_find) THEN
         xhelp1 = tof12(1,tof12_i,itdc)+tof12(2,tof12_i,itdc)
         xhelp2 = tof21(1,tof21_i,itdc)+tof21(2,tof21_i,itdc)
         ds = xhelp1-xhelp2
         ihelp=(tof12_i-1)*2+tof21_i
         c1 = k_S12S21(1,ihelp)
         c2 = k_S12S21(2,ihelp)
	 if ((xhelp1.lt.8000.).and.(xhelp2.lt.8000)) 
     &   beta_a(11) = c2/(cos(theta12)*(ds-c1))
      ENDIF

C     S12 - S22
      IF (tof12_i.GT.none_find.AND.tof22_i.GT.none_find) THEN
         xhelp1 = tof12(1,tof12_i,itdc)+tof12(2,tof12_i,itdc)
         xhelp2 = tof22(1,tof22_i,itdc)+tof22(2,tof22_i,itdc)
         ds = xhelp1-xhelp2
         ihelp=(tof12_i-1)*2+tof22_i
         c1 = k_S12S22(1,ihelp)
         c2 = k_S12S22(2,ihelp)
	 if ((xhelp1.lt.8000.).and.(xhelp2.lt.8000)) 
     &   beta_a(12) = c2/(cos(theta12)*(ds-c1))
      ENDIF
      
C-------

      icount=0
      sw=0.
      sxw=0.
      beta_mean=100.

      do i=1,12
         if ((beta_a(i).gt.-1.5).and.(beta_a(i).lt.1.5)) then
            icount= icount+1
            if (i.le.4) w_i=1./(0.13**2.)
            if ((i.ge.5).and.(i.le.8)) w_i=1./(0.16**2.)
	    if (i.ge.9) w_i=1./(0.25**2.)     ! to be checked
            sxw=sxw + beta_a(i)*w_i
            sw =sw + w_i
         endif
      enddo
      
      if (icount.gt.0) beta_mean=sxw/sw
      beta_a(13) = beta_mean

      RETURN
      END


