ToFLevel2/src/tofl2com.for


C******************************************************************************
C
C  08-12-06 WM: adc_c-bug :  The raw ADc value was multiplied with cos(theta)
C  and AFTER that there was an if statement "if tof32(right,i,iadc) < 4095"
C
C  jan-07 GF: ADCflags(4,12) inserted to flag artificial ADC values
C  jan-07 WM: artificial ADC values created using attenuation calibration
C  jan-07 WM: modified xtofpos flag "101". xtofpos must be inside physical
C             dimension of the paddle +/- 10 cm
C  jan-07 WM: if xtofpos=101 then this paddle is not used for beta
C             calculation
C  jan-07 WM: the definition for a "hit" is changed: Now we must have a
C             valid TDC signal on both sides
C  jan-07 WM: flag for PMTs #10 and #35 added, TDC=819 due to bit-shift
C  jan-07 WM: bug removed: in some cases tdc_tw was calculated due to a
C             leftover "xhelp" value
C  apr-07 WM: attenuation fit curve is now a double exponential fit
C             conversion from raw ADC to pC using calibration function
C             variables xtr_tof and ytr_tof inserted (filled with default) 
C  jan-08 WM: Major Update: Time Walk correction introduced
C             Additionalyl we use the information from the "check_charge"
C             function to fill artificial ADC values and make small corrections
C             to the k1-parameter (for Z>2)
C  feb-08 WM: Calculation of beta(13) changed: First a mean beta is calculated,
C             then in a second step we check the residuals of the single 
C             measurements, reject if > 10 sigma, calculate chi2 and "quality"
C             beta is taken as good if chi2<20 and quality>10
C  mar-08 WM: Call to "newbeta" changed, now a flag tells the function if the
C             call comes from "tofl2com" or form "toftrack"
C  mar-08 WM: Bug found in dEdx if check_charge>1
C  oct-08 WM: Calculation of zenith angle debugged, sometimes strange values
C             were possible
C  nov-09 WM: the dEdx part ("adctof_c") moved to the new dEdx routine from Napoli
C  feb-10 WM: k1 values now for Z=1, Z=2, Z>2, k2 values are fix
C******************************************************************************

      INTEGER FUNCTION TOFL2COM()
c     
      IMPLICIT NONE
C     
      include  'input_tof.txt'
      include  'output_tof.txt'
      include  'tofcomm.txt'

      INTEGER icounter
      DATA icounter / 0/

      LOGICAL check
      REAL secure

      INTEGER j,hitvec(6)

      REAL dx,dy,dr,ds
      REAL yhelp,yhelp1,yhelp2,xhelp,xhelp1,xhelp2
      REAL c1,c2
      REAL dist

C      REAL sw,sxw,w_i
C      INTEGER icount
C      REAL beta_mean

      INTEGER tof11_j,tof21_j,tof31_j
      INTEGER tof12_j,tof22_j,tof32_j

c     value for status of each PM-data
c     first index  : 1 = left, 2 = right
c     second index : 1... number of paddle
      INTEGER tof11_event(2,8),tof12_event(2,6)
      INTEGER tof21_event(2,2),tof22_event(2,2)
      INTEGER tof31_event(2,3),tof32_event(2,3)


      REAL y_coor_lin11c(8,2),x_coor_lin12c(6,2)
      REAL x_coor_lin21c(2,2),y_coor_lin22c(2,2)
      REAL y_coor_lin31c(3,2),x_coor_lin32c(3,2)

      DATA y_coor_lin11c(1,1),y_coor_lin11c(1,2) /-20.66,-2.497/
      DATA y_coor_lin11c(2,1),y_coor_lin11c(2,2) /-9.10, -2.52/
      DATA y_coor_lin11c(3,1),y_coor_lin11c(3,2) /-24.07,-2.12/
      DATA y_coor_lin11c(4,1),y_coor_lin11c(4,2) /-13.40,-2.47/
      DATA y_coor_lin11c(5,1),y_coor_lin11c(5,2) /-31.07,-2.32/
      DATA y_coor_lin11c(6,1),y_coor_lin11c(6,2) /-21.69,-2.63/
      DATA y_coor_lin11c(7,1),y_coor_lin11c(7,2) /-12.37,-2.65/
      DATA y_coor_lin11c(8,1),y_coor_lin11c(8,2) /-10.81,-3.15/

      DATA x_coor_lin12c(1,1),x_coor_lin12c(1,2) /12.96, -2.65/
      DATA x_coor_lin12c(2,1),x_coor_lin12c(2,2) /17.12,-2.44/
      DATA x_coor_lin12c(3,1),x_coor_lin12c(3,2) /7.26, -1.98/
      DATA x_coor_lin12c(4,1),x_coor_lin12c(4,2) /-22.52,-2.27/
      DATA x_coor_lin12c(5,1),x_coor_lin12c(5,2) /-18.54,-2.28/
      DATA x_coor_lin12c(6,1),x_coor_lin12c(6,2) /-7.67,-2.15/

      DATA x_coor_lin21c(1,1),x_coor_lin21c(1,2) /22.56,-1.56/
      DATA x_coor_lin21c(2,1),x_coor_lin21c(2,2) /13.94,-1.56/

      DATA y_coor_lin22c(1,1),y_coor_lin22c(1,2) /-24.24,-2.23/
      DATA y_coor_lin22c(2,1),y_coor_lin22c(2,2) /-45.99,-1.68/

      DATA y_coor_lin31c(1,1),y_coor_lin31c(1,2) /-22.99,-3.54/
      DATA y_coor_lin31c(2,1),y_coor_lin31c(2,2) /-42.28,-4.10/
      DATA y_coor_lin31c(3,1),y_coor_lin31c(3,2) /-41.29,-3.69/

      DATA x_coor_lin32c(1,1),x_coor_lin32c(1,2) /0.961, -3.22/
      DATA x_coor_lin32c(2,1),x_coor_lin32c(2,2) /4.98,-3.48/
      DATA x_coor_lin32c(3,1),x_coor_lin32c(3,2) /-22.08,-3.37/

      
      REAL theta13

c      DOUBLE PRECISION ZTOF(6)
      REAL ZTOF(6) !EM GCC4.7
      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

      REAL hepratio      

      INTEGER ihelp
      REAL xkorr,btemp(12)

      REAL atten,pc_adc,check_charge,newbeta

      INTEGER IZ


      INTEGER ifst 
      DATA  ifst /0/

C---------------------------------------
C     
C     Begin !
C     
      TOFL2COM = 0
C     
C     CALCULATE COMMON VARIABLES
C     
C-------------------------------------------------------------------

        if (ifst.eq.0) then

        ifst=1

C     amplitude has to be 'secure' higher than pedestal for an adc event
       secure = 2.

C     ratio between helium and proton ca. 4
       hepratio = 4.  ! 
       offset = 1
       slope = 2
       left = 1
       right = 2
       none_ev = 0
       none_find = 0
       tdc_ev = 1
       adc_ev = 1
       itdc = 1
       iadc = 2

        ENDIF
C---------------------------------------------------------------------

      icounter = icounter + 1


      do i=1,13
         betatof_a(i) = 100.          ! As in "troftrk.for"
      enddo

      do i=1,6
         hitvec(i) = -1
      enddo

      do i=1,4
         do j=1,12
            adctof_c(i,j) = 1000.
         enddo
      enddo


      do i=1,4
         do j=1,12
            tdc_c(i,j) = 4095.
         enddo
      enddo


      do i=1,12
         do j=1,4
            tofmask(j,i) = 0
         enddo
      enddo


c gf adc falg:
      do i=1,4
         do j=1,12
            adcflagtof(i,j) = 0
         enddo
      enddo

c gf tdc falg:
      do i=1,4
         do j=1,12
            tdcflagtof(i,j) = 0
         enddo
      enddo


C---  Fill xtr_tof and ytr_tof: positions from tracker at ToF layers
C---  since this is standalone ToF fill with default values
      do j=1,6
      xtr_tof(j) = 101.
      ytr_tof(j) = 101.
      enddo

c the calibration files are read in the main program from xxx_tofcalib.rz

c-------------------------get ToF data --------------------------------

c     put the adc and tdc values from ntuple into tofxx(i,j,k) variables
c     adc valueas are then pC

      do j=1,8
         tof11(1,j,2) = pc_adc(adc(ch11a(j),hb11a(j)))
         tof11(2,j,2) = pc_adc(adc(ch11b(j),hb11b(j)))
         tof11(1,j,1) = (tdc(ch11a(j),hb11a(j)))
         tof11(2,j,1) = (tdc(ch11b(j),hb11b(j)))
      enddo


      do j=1,6
         tof12(1,j,2) = pc_adc(adc(ch12a(j),hb12a(j)))
         tof12(2,j,2) = pc_adc(adc(ch12b(j),hb12b(j)))
         tof12(1,j,1) = (tdc(ch12a(j),hb12a(j)))
         tof12(2,j,1) = (tdc(ch12b(j),hb12b(j)))
      enddo

      do j=1,2
         tof21(1,j,2) = pc_adc(adc(ch21a(j),hb21a(j)))
         tof21(2,j,2) = pc_adc(adc(ch21b(j),hb21b(j)))
         tof21(1,j,1) = (tdc(ch21a(j),hb21a(j)))
         tof21(2,j,1) = (tdc(ch21b(j),hb21b(j)))
      enddo

      do j=1,2
         tof22(1,j,2) = pc_adc(adc(ch22a(j),hb22a(j)))
         tof22(2,j,2) = pc_adc(adc(ch22b(j),hb22b(j)))
         tof22(1,j,1) = (tdc(ch22a(j),hb22a(j)))
         tof22(2,j,1) = (tdc(ch22b(j),hb22b(j)))
      enddo

      do j=1,3
         tof31(1,j,2) = pc_adc(adc(ch31a(j),hb31a(j)))
         tof31(2,j,2) = pc_adc(adc(ch31b(j),hb31b(j)))
         tof31(1,j,1) = (tdc(ch31a(j),hb31a(j)))
         tof31(2,j,1) = (tdc(ch31b(j),hb31b(j)))
      enddo

      do j=1,3
         tof32(1,j,2) = pc_adc(adc(ch32a(j),hb32a(j)))
         tof32(2,j,2) = pc_adc(adc(ch32b(j),hb32b(j)))
         tof32(1,j,1) = (tdc(ch32a(j),hb32a(j)))
         tof32(2,j,1) = (tdc(ch32b(j),hb32b(j)))
      enddo

C----------------------------------------------------------------------

      DO i = 1,8
         if (abs(tof11(1,i,itdc)).gt.10000.) tof11(1,i,itdc)= 10000.
         if (abs(tof11(2,i,itdc)).gt.10000.) tof11(2,i,itdc)= 10000.
         if (abs(tof11(1,i,iadc)).gt.10000.) tof11(1,i,iadc)= 10000.
         if (abs(tof11(2,i,iadc)).gt.10000.) tof11(2,i,iadc)= 10000.
      ENDDO

      DO i = 1,6
         if (abs(tof12(1,i,itdc)).gt.10000.) tof12(1,i,itdc)= 10000.
         if (abs(tof12(2,i,itdc)).gt.10000.) tof12(2,i,itdc)= 10000.
         if (abs(tof12(1,i,iadc)).gt.10000.) tof12(1,i,iadc)= 10000.
         if (abs(tof12(2,i,iadc)).gt.10000.) tof12(2,i,iadc)= 10000.
      ENDDO


      DO i = 1,2
         if (abs(tof21(1,i,itdc)).gt.10000.) tof21(1,i,itdc)= 10000.
         if (abs(tof21(2,i,itdc)).gt.10000.) tof21(2,i,itdc)= 10000.
         if (abs(tof21(1,i,iadc)).gt.10000.) tof21(1,i,iadc)= 10000.
         if (abs(tof21(2,i,iadc)).gt.10000.) tof21(2,i,iadc)= 10000.
      ENDDO

      DO i = 1,2
         if (abs(tof22(1,i,itdc)).gt.10000.) tof22(1,i,itdc)= 10000.
         if (abs(tof22(2,i,itdc)).gt.10000.) tof22(2,i,itdc)= 10000.
         if (abs(tof22(1,i,iadc)).gt.10000.) tof22(1,i,iadc)= 10000.
         if (abs(tof22(2,i,iadc)).gt.10000.) tof22(2,i,iadc)= 10000.
      ENDDO

      DO i = 1,3
         if (abs(tof31(1,i,itdc)).gt.10000.) tof31(1,i,itdc)= 10000.
         if (abs(tof31(2,i,itdc)).gt.10000.) tof31(2,i,itdc)= 10000.
         if (abs(tof31(1,i,iadc)).gt.10000.) tof31(1,i,iadc)= 10000.
         if (abs(tof31(2,i,iadc)).gt.10000.) tof31(2,i,iadc)= 10000.
      ENDDO

      DO i = 1,3
         if (abs(tof32(1,i,itdc)).gt.10000.) tof32(1,i,itdc)= 10000.
         if (abs(tof32(2,i,itdc)).gt.10000.) tof32(2,i,itdc)= 10000.
         if (abs(tof32(1,i,iadc)).gt.10000.) tof32(1,i,iadc)= 10000.
         if (abs(tof32(2,i,iadc)).gt.10000.) tof32(2,i,iadc)= 10000.
      ENDDO

C----------------------------------------------------------------------
C------------------  set ADC & TDC flag = 0    ------------------------
C----------------------------------------------------------------------

      do j=1,8
      if (adc(ch11a(j),hb11a(j)).LT.4096)adcflagtof(ch11a(j),hb11a(j))=0
      if (adc(ch11b(j),hb11b(j)).LT.4096)adcflagtof(ch11b(j),hb11b(j))=0
      if (tdc(ch11a(j),hb11a(j)).LT.4096)tdcflagtof(ch11a(j),hb11a(j))=0
      if (tdc(ch11b(j),hb11b(j)).LT.4096)tdcflagtof(ch11b(j),hb11b(j))=0
      enddo
      do j=1,6
      if (adc(ch12a(j),hb12a(j)).LT.4096)adcflagtof(ch12a(j),hb12a(j))=0
      if (adc(ch12b(j),hb12b(j)).LT.4096)adcflagtof(ch12b(j),hb12b(j))=0
      if (tdc(ch12a(j),hb12a(j)).LT.4096)tdcflagtof(ch12a(j),hb12a(j))=0
      if (tdc(ch12b(j),hb12b(j)).LT.4096)tdcflagtof(ch12b(j),hb12b(j))=0
      enddo
      do j=1,2
      if (adc(ch21a(j),hb21a(j)).LT.4096)adcflagtof(ch21a(j),hb21a(j))=0
      if (adc(ch21b(j),hb21b(j)).LT.4096)adcflagtof(ch21b(j),hb21b(j))=0
      if (tdc(ch21a(j),hb21a(j)).LT.4096)tdcflagtof(ch21a(j),hb21a(j))=0
      if (tdc(ch21b(j),hb21b(j)).LT.4096)tdcflagtof(ch21b(j),hb21b(j))=0
      enddo
      do j=1,2
      if (adc(ch22a(j),hb22a(j)).LT.4096)adcflagtof(ch22a(j),hb22a(j))=0
      if (adc(ch22b(j),hb22b(j)).LT.4096)adcflagtof(ch22b(j),hb22b(j))=0
      if (tdc(ch22a(j),hb22a(j)).LT.4096)tdcflagtof(ch22a(j),hb22a(j))=0
      if (tdc(ch22b(j),hb22b(j)).LT.4096)tdcflagtof(ch22b(j),hb22b(j))=0
      enddo
      do j=1,3
      if (adc(ch31a(j),hb31a(j)).LT.4096)adcflagtof(ch31a(j),hb31a(j))=0
      if (adc(ch31b(j),hb31b(j)).LT.4096)adcflagtof(ch31b(j),hb31b(j))=0
      if (tdc(ch31a(j),hb31a(j)).LT.4096)tdcflagtof(ch31a(j),hb31a(j))=0
      if (tdc(ch31b(j),hb31b(j)).LT.4096)tdcflagtof(ch31b(j),hb31b(j))=0
      enddo
      do j=1,3
      if (adc(ch32a(j),hb32a(j)).LT.4096)adcflagtof(ch32a(j),hb32a(j))=0
      if (adc(ch32b(j),hb32b(j)).LT.4096)adcflagtof(ch32b(j),hb32b(j))=0
      if (tdc(ch32a(j),hb32a(j)).LT.4096)tdcflagtof(ch32a(j),hb32a(j))=0
      if (tdc(ch32b(j),hb32b(j)).LT.4096)tdcflagtof(ch32b(j),hb32b(j))=0
      enddo

C----------------------------------------------------------------
C---------- Check PMTs 10 and 35 for strange TDC values----------
C----------------------------------------------------------------

C---- S116A TDC=819
       if (tof11(1,6,1).EQ.819) then
             tof11(1,6,1) = 4095
             tdcflagtof(ch11a(6),hb11a(6))=2
       endif
   
C---- S222B TDC=819
       if (tof22(2,2,1).EQ.819) then
             tof22(2,2,1) = 4095
             tdcflagtof(ch22b(2),hb22b(2))=2
       endif
 
C----------------------------------------------------------------
C------------   Check Paddles for hits    -----------------------
C------  a "hit" means TDC values<4095 on both sides ------------
C----------------------------------------------------------------

C     upper tof  S11
      DO i = 1,8

         DO j = 1,2
            tof11_event(j,i) = none_ev
            IF ((tof11(j,i,itdc).LT.2000).AND.(tof11(j,i,itdc).GT.100))
     +           tof11_event(j,i) = tof11_event(j,i) + tdc_ev
         ENDDO
      ENDDO

c     find single paddle in upper tof with tdc and adc signal
      tof11_i = none_find
      tof11_j = none_find
      check = .TRUE.
      DO i = 1, 8
         IF ((tof11_event(left,i).GE.1).AND.(tof11_event(right,i).GE.1))
     +        THEN
c     check if an other paddle has also an event - then set flag
            tof11_j = tof11_j + 2**(i-1)
            IF (check.EQV..TRUE.) THEN
               IF (tof11_i.EQ.none_find) THEN
                  tof11_i = i
               ELSE
                  tof11_i = -1
                  check = .FALSE.
               ENDIF
            ENDIF
         ENDIF
      ENDDO


C     upper tof  S12
      DO i = 1,6
         DO j = 1,2
            tof12_event(j,i) = none_ev
            IF ((tof12(j,i,itdc).LT.2000).AND.(tof12(j,i,itdc).GT.100))
     +           tof12_event(j,i) = tof12_event(j,i) + tdc_ev
         ENDDO
      ENDDO

c     find single paddle in upper tof with tdc and adc signal
      tof12_i = none_find
      tof12_j = none_find
      check = .TRUE.
      DO i = 1, 6
         IF ((tof12_event(left,i).GE.1).AND.(tof12_event(right,i).GE.1))
     +        THEN
c     check if an other paddle has also an event - then set flag
            tof12_j = tof12_j + 2**(i-1)
            IF (check.EQV..TRUE.) THEN
               IF (tof12_i.EQ.none_find) THEN
                  tof12_i = i
               ELSE
                  tof12_i = -1
                  check = .FALSE.
               ENDIF
            ENDIF
         ENDIF
      ENDDO


C     middle tof  S21
      DO i = 1,2
         DO j = 1,2
            tof21_event(j,i) = none_ev
            IF ((tof21(j,i,itdc).LT.2000).AND.(tof21(j,i,itdc).GT.100))
     +           tof21_event(j,i) = tof21_event(j,i) + tdc_ev
         ENDDO
      ENDDO

c     find single paddle in upper tof with tdc and adc signal
      tof21_i = none_find
      tof21_j = none_find
      check = .TRUE.
      DO i = 1, 2
         IF ((tof21_event(left,i).GE.1).AND.(tof21_event(right,i).GE.1))
     +        THEN
c     check if an other paddle has also an event - then set flag
            tof21_j = tof21_j + 2**(i-1)
            IF (check.EQV..TRUE.) THEN
               IF (tof21_i.EQ.none_find) THEN
                  tof21_i = i
               ELSE
                  tof21_i = -1
                  check = .FALSE.
               ENDIF
            ENDIF
         ENDIF
      ENDDO

C     middle tof  S22
      DO i = 1,2
         DO j = 1,2
            tof22_event(j,i) = none_ev
            IF ((tof22(j,i,itdc).LT.2000).AND.(tof22(j,i,itdc).GT.100))
     +           tof22_event(j,i) = tof22_event(j,i) + tdc_ev
         ENDDO
      ENDDO

c     find single paddle in upper tof with tdc and adc signal
      tof22_i = none_find
      tof22_j = none_find
      check = .TRUE.
      DO i = 1, 2
         IF ((tof22_event(left,i).GE.1).AND.(tof22_event(right,i).GE.1))
     +        THEN
c     check if an other paddle has also an event - then set flag
            tof22_j = tof22_j + 2**(i-1)
            IF (check.EQV..TRUE.) THEN
               IF (tof22_i.EQ.none_find) THEN
                  tof22_i = i
               ELSE
                  tof22_i = -1
                  check = .FALSE.
               ENDIF
            ENDIF
         ENDIF
      ENDDO


C     bottom tof  S31
      DO i = 1,3
         DO j = 1,2
            tof31_event(j,i) = none_ev
            IF ((tof31(j,i,itdc).LT.2000).AND.(tof31(j,i,itdc).GT.100))
     +           tof31_event(j,i) = tof31_event(j,i) + tdc_ev
         ENDDO
      ENDDO

c     find single paddle in upper tof with tdc and adc signal
      tof31_i = none_find
      tof31_j = none_find
      check = .TRUE.
      DO i = 1, 3
         IF ((tof31_event(left,i).GE.1).AND.(tof31_event(right,i).GE.1))
     +        THEN
c     check if an other paddle has also an event - then set flag
            tof31_j = tof31_j + 2**(i-1)
            IF (check.EQV..TRUE.) THEN
               IF (tof31_i.EQ.none_find) THEN
                  tof31_i = i
               ELSE
                  tof31_i = -1
                  check = .FALSE.
               ENDIF
            ENDIF
         ENDIF
      ENDDO

C     bottom tof  S32
      DO i = 1,3
         DO j = 1,2
            tof32_event(j,i) = none_ev
            IF ((tof32(j,i,itdc).LT.2000).AND.(tof32(j,i,itdc).GT.100))
     +           tof32_event(j,i) = tof32_event(j,i) + tdc_ev
         ENDDO
      ENDDO

c     find single paddle in upper tof with tdc and adc signal
      tof32_i = none_find
      tof32_j = none_find
      check = .TRUE.
      DO i = 1, 3
         IF ((tof32_event(left,i).GE.1).AND.(tof32_event(right,i).GE.1))
     +        THEN
c     check if an other paddle has also an event - then set flag
            tof32_j = tof32_j + 2**(i-1)
            IF (check.EQV..TRUE.) THEN
               IF (tof32_i.EQ.none_find) THEN
                  tof32_i = i
               ELSE
                  tof32_i = -1
                  check = .FALSE.
               ENDIF
            ENDIF
         ENDIF
       ENDDO

       do i=1,6
         tof_i_flag(i)=0
         tof_j_flag(i)=0
       enddo 
 
       tof_i_flag(1)=tof11_i
       tof_i_flag(2)=tof12_i
       tof_i_flag(3)=tof21_i
       tof_i_flag(4)=tof22_i
       tof_i_flag(5)=tof31_i
       tof_i_flag(6)=tof32_i

       tof_j_flag(1)=tof11_j
       tof_j_flag(2)=tof12_j
       tof_j_flag(3)=tof21_j 
       tof_j_flag(4)=tof22_j
       tof_j_flag(5)=tof31_j
       tof_j_flag(6)=tof32_j

       hitvec(1)=tof11_i
       hitvec(2)=tof12_i
       hitvec(3)=tof21_i
       hitvec(4)=tof22_i
       hitvec(5)=tof31_i
       hitvec(6)=tof32_i


C------------------------------------------------------------------
C--  calculate track position in paddle using timing difference
C--  this calculation is preliminary and uses some standard 
C--  calibration values, but we need to find a rough position to
C--  be able to calculate artificial ADC values (needed for the 
C--  timewalk...
C------------------------------------------------------------------

       do i=1,3
         xtofpos(i)=100.
         ytofpos(i)=100.
       enddo 

C-----------------------------S1 --------------------------------

      IF (tof11_i.GT.none_find) THEN
         ytofpos(1)  = ((tof11(1,tof11_i,itdc)-tof11(2,tof11_i,itdc))/2.
     +   -y_coor_lin11c(tof11_i,offset))/y_coor_lin11c(tof11_i,slope)
      endif

      IF (tof12_i.GT.none_find) THEN
         xtofpos(1)  = ((tof12(1,tof12_i,itdc)-tof12(2,tof12_i,itdc))/2.
     +   -x_coor_lin12c(tof12_i,offset))/x_coor_lin12c(tof12_i,slope)
      endif


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

      IF (tof21_i.GT.none_find) THEN
         xtofpos(2) = ((tof21(1,tof21_i,itdc)-tof21(2,tof21_i,itdc))/2.
     +    -x_coor_lin21c(tof21_i,offset))/x_coor_lin21c(tof21_i,slope)
      endif

      IF (tof22_i.GT.none_find) THEN
         ytofpos(2) = ((tof22(1,tof22_i,itdc)-tof22(2,tof22_i,itdc))/2.
     +    -y_coor_lin22c(tof22_i,offset))/y_coor_lin22c(tof22_i,slope)
      endif
      

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

      IF (tof31_i.GT.none_find) THEN
         ytofpos(3)  = ((tof31(1,tof31_i,itdc)-tof31(2,tof31_i,itdc))/2.
     +    -y_coor_lin31c(tof31_i,offset))/y_coor_lin31c(tof31_i,slope)
      endif

      IF (tof32_i.GT.none_find) THEN
         xtofpos(3)  = ((tof32(1,tof32_i,itdc)-tof32(2,tof32_i,itdc))/2.
     +    -x_coor_lin32c(tof32_i,offset))/x_coor_lin32c(tof32_i,slope)
      endif


C----------------------------------------------------------------------
C---------------------  zenith angle theta  ---------------------------
C----------------------------------------------------------------------

       xhelp1=0.
       if (tof11_i.GT.none_find) xhelp1=tof11_x(tof11_i)
       if (xtofpos(1).lt.100)  xhelp1=xtofpos(1)

       yhelp1=0.
       if (tof12_i.GT.none_find) yhelp1=tof12_y(tof12_i)
       if (ytofpos(1).lt.100)  yhelp1=ytofpos(1)


       yhelp2=0.
       if (tof32_i.GT.none_find) yhelp2=tof32_y(tof32_i)
       if (ytofpos(3).lt.100)  yhelp2=ytofpos(3)

       xhelp2=0.
       if (tof31_i.GT.none_find) xhelp2=tof31_x(tof31_i)
       if (xtofpos(3).lt.100)  xhelp2=xtofpos(3)


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

       dx  = xhelp1 - xhelp2
       dy  = yhelp1 - yhelp2
       dr = sqrt(dx*dx+dy*dy)
       theta13 = atan(dr/tofarm13)


C----------------------------------------------------------------------
C--- check charge:
C--- if Z=2 we should use the attenuation curve for helium to 
C--- fill the artificail ADC values and NOT divide by "hepratio"
C--- if Z>2 we should do a correction to 
C--- the k1 constants in the beta calculation
C----------------------------------------------------------------------

         iz = int(check_charge(theta13,hitvec))
C         write(*,*) 'charge in tofl2com',iz

C--------------------------------------------------------------------
C---- if TDCleft.and.TDCright and NO ADC insert artificial ADC 
C---- values 
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/     ! paddles in different order
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----------------------------  S1 -------------------------------------

c       yhelp=0.
       yhelp=100.  ! WM
       if (tof12_i.GT.none_find) yhelp=tof12_y(tof12_i)
       if (ytofpos(1).lt.100)  yhelp=ytofpos(1)

       IF (tof11_i.GT.none_find.AND.abs(yhelp).lt.100) THEN
         i = tof11_i
         if (adc(ch11a(i),hb11a(i)).eq.4095) then
            xkorr = atten(left,11,i,yhelp)
            if (iz.le.1) xkorr=xkorr/hepratio
            tof11(left,i,iadc)=xkorr/cos(theta13)
            adcflagtof(ch11a(i),hb11a(i)) = 1
         endif
         if (adc(ch11b(i),hb11b(i)).eq.4095) then
            xkorr = atten(right,11,i,yhelp)
            if (iz.le.1) xkorr=xkorr/hepratio
            tof11(right,i,iadc)=xkorr/cos(theta13)
            adcflagtof(ch11b(i),hb11b(i)) = 1
         endif
       ENDIF

c       xhelp=0.
       xhelp=100.  ! WM
       if (tof11_i.GT.none_find) xhelp=tof11_x(tof11_i)
       if (xtofpos(1).lt.100)  xhelp=xtofpos(1)

       IF (tof12_i.GT.none_find.AND.abs(xhelp).lt.100) THEN
         i = tof12_i
         if (adc(ch12a(i),hb12a(i)).eq.4095) then
            xkorr = atten(left,12,i,xhelp)
            if (iz.le.1) xkorr=xkorr/hepratio
            tof12(left,i,iadc) = xkorr/cos(theta13)
            adcflagtof(ch12a(i),hb12a(i)) = 1
         endif
         if (adc(ch12b(i),hb12b(i)).eq.4095) then
            xkorr = atten(right,12,i,xhelp)
            if (iz.le.1) xkorr=xkorr/hepratio
            tof12(right,i,iadc) = xkorr/cos(theta13)
            adcflagtof(ch12b(i),hb12b(i)) = 1
         endif
       ENDIF

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

c       xhelp=0.
       xhelp=100.   ! WM
       if (tof22_i.GT.none_find) xhelp=tof22_x(tof22_i)
       if (xtofpos(2).lt.100)  xhelp=xtofpos(2)

       IF (tof21_i.GT.none_find.AND.abs(xhelp).lt.100) THEN
         i = tof21_i
         if (adc(ch21a(i),hb21a(i)).eq.4095) then
            xkorr = atten(left,21,i,xhelp)
            if (iz.le.1) xkorr=xkorr/hepratio
            tof21(left,i,iadc) = xkorr/cos(theta13)
            adcflagtof(ch21a(i),hb21a(i)) = 1
         endif
         if (adc(ch21b(i),hb21b(i)).eq.4095) then
            xkorr = atten(right,21,i,xhelp)
            if (iz.le.1) xkorr=xkorr/hepratio
            tof21(right,i,iadc) = xkorr/cos(theta13)
            adcflagtof(ch21b(i),hb21b(i)) = 1
         endif
       ENDIF


c       yhelp=0.
       yhelp=100.   ! WM
       if (tof21_i.GT.none_find) yhelp=tof21_y(tof21_i)
       if (ytofpos(2).lt.100)  yhelp=ytofpos(2)

       IF (tof22_i.GT.none_find.AND.abs(yhelp).lt.100) THEN
         i = tof22_i
         if (adc(ch22a(i),hb22a(i)).eq.4095) then
            xkorr = atten(left,22,i,yhelp)
            if (iz.le.1) xkorr=xkorr/hepratio
            tof22(left,i,iadc) = xkorr/cos(theta13)
            adcflagtof(ch22a(i),hb22a(i)) = 1
         endif
         if (adc(ch22b(i),hb22b(i)).eq.4095) then
            xkorr = atten(right,22,i,yhelp)
            if (iz.le.1) xkorr=xkorr/hepratio
            tof22(right,i,iadc) = xkorr/cos(theta13)
            adcflagtof(ch22b(i),hb22b(i)) = 1
         endif
       ENDIF

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

c       yhelp=0.
       yhelp=100.  ! WM
       if (tof32_i.GT.none_find) yhelp=tof32_y(tof32_i)
       if (ytofpos(3).lt.100)  yhelp=ytofpos(3)

       IF (tof31_i.GT.none_find.AND.abs(yhelp).lt.100) THEN
         i = tof31_i
         if (adc(ch31a(i),hb31a(i)).eq.4095) then
            xkorr = atten(left,31,i,yhelp)
            if (iz.le.1) xkorr=xkorr/hepratio
            tof31(left,i,iadc) = xkorr/cos(theta13)
            adcflagtof(ch31a(i),hb31a(i)) = 1
         endif
         if (adc(ch31b(i),hb31b(i)).eq.4095) then
            xkorr = atten(right,31,i,yhelp)
            if (iz.le.1) xkorr=xkorr/hepratio
            tof31(right,i,iadc) = xkorr/cos(theta13)
            adcflagtof(ch31b(i),hb31b(i)) = 1
         endif
       ENDIF

c       xhelp=0.
       xhelp=100.   ! WM
       if (tof31_i.GT.none_find) xhelp=tof31_x(tof31_i)
       if (xtofpos(3).lt.100)  xhelp=xtofpos(3)

       IF (tof32_i.GT.none_find.AND.abs(xhelp).lt.100) THEN
         i = tof32_i
         if (adc(ch32a(i),hb32a(i)).eq.4095) then
            xkorr = atten(left,32,i,xhelp)
            if (iz.le.1) xkorr=xkorr/hepratio
            tof32(left,i,iadc) = xkorr/cos(theta13)
            adcflagtof(ch32a(i),hb32a(i)) = 1
         endif
         if (adc(ch32b(i),hb32b(i)).eq.4095) then
            xkorr = atten(right,32,i,xhelp)
            if (iz.le.1) xkorr=xkorr/hepratio
            tof32(right,i,iadc) = xkorr/cos(theta13)
            adcflagtof(ch32b(i),hb32b(i)) = 1
         endif
       ENDIF


C-------------------------------------------------------------------
C--------------------Time walk correction  -------------------------
C-------------------------------------------------------------------
C-------------------------------------------------------------------
C Now there is for each hitted paddle a TDC and ADC value, if the
C TDC was < 4095.
C There might be also TDC-ADC pairs in paddles not hitted

C-------------------------------------------------------------------
C If we have multiple paddles hit, so that no artificial ADC value
C is created, we set the raw TDC value as "tdc_c"
C-------------------------------------------------------------------
c 
c       do i=1,4
c         do j=1,12
c            tdc_c(i,j) = tdc(i,j)
c         enddo
c       enddo
c
C----  Let's correct the raw TDC value with the time walk  ---------

      DO i=1,8
         if ((tdc(ch11a(i),hb11a(i)).lt.4095).and.
     &             (tof11(left,i,iadc).lt.3786)) THEN
         xhelp = tw11(left,i)/(tof11(left,i,iadc)**0.5)
         tof11(left,i,itdc) = tof11(left,i,itdc) + xhelp
         tdc_c(ch11a(i),hb11a(i))=tof11(left,i,itdc)
                                              ENDIF

         if ((tdc(ch11b(i),hb11b(i)).lt.4095).and.
     &             (tof11(right,i,iadc).lt.3786)) THEN
         xhelp = tw11(right,i)/(tof11(right,i,iadc)**0.5)
         tof11(right,i,itdc) = tof11(right,i,itdc) + xhelp
         tdc_c(ch11b(i),hb11b(i))=tof11(right,i,itdc)
                                             ENDIF
      ENDDO


      DO i=1,6
         if ((tdc(ch12a(i),hb12a(i)).lt.4095).and.
     &             (tof12(left,i,iadc).lt.3786)) THEN
         xhelp = tw12(left,i)/(tof12(left,i,iadc)**0.5)
         tof12(left,i,itdc) = tof12(left,i,itdc) + xhelp
         tdc_c(ch12a(i),hb12a(i))=tof12(left,i,itdc)
                                              ENDIF

         if ((tdc(ch12b(i),hb12b(i)).lt.4095).and.
     &             (tof12(right,i,iadc).lt.3786)) THEN
         xhelp = tw12(right,i)/(tof12(right,i,iadc)**0.5)
         tof12(right,i,itdc) = tof12(right,i,itdc) + xhelp
         tdc_c(ch12b(i),hb12b(i))=tof12(right,i,itdc)
                                             ENDIF
      ENDDO

C----
      DO I=1,2
         if ((tdc(ch21a(i),hb21a(i)).lt.4095).and.
     &             (tof21(left,i,iadc).lt.3786)) THEN
         xhelp = tw21(left,i)/(tof21(left,i,iadc)**0.5)
         tof21(left,i,itdc) = tof21(left,i,itdc) + xhelp
         tdc_c(ch21a(i),hb21a(i))=tof21(left,i,itdc)
                                              ENDIF

         if ((tdc(ch21b(i),hb21b(i)).lt.4095).and.
     &             (tof21(right,i,iadc).lt.3786)) THEN
         xhelp = tw21(right,i)/(tof21(right,i,iadc)**0.5)
         tof21(right,i,itdc) = tof21(right,i,itdc) + xhelp
         tdc_c(ch21b(i),hb21b(i))=tof21(right,i,itdc)
                                             ENDIF
      ENDDO

      DO I=1,2
         if ((tdc(ch22a(i),hb22a(i)).lt.4095).and.
     &             (tof22(left,i,iadc).lt.3786)) THEN
         xhelp = tw22(left,i)/(tof22(left,i,iadc)**0.5)
         tof22(left,i,itdc) = tof22(left,i,itdc) + xhelp
         tdc_c(ch22a(i),hb22a(i))=tof22(left,i,itdc)
                                              ENDIF

         if ((tdc(ch22b(i),hb22b(i)).lt.4095).and.
     &             (tof22(right,i,iadc).lt.3786)) THEN
         xhelp = tw22(right,i)/(tof22(right,i,iadc)**0.5)
         tof22(right,i,itdc) = tof22(right,i,itdc) + xhelp
         tdc_c(ch22b(i),hb22b(i))=tof22(right,i,itdc)
                                             ENDIF
      ENDDO

C----
      DO I=1,3
         if ((tdc(ch31a(i),hb31a(i)).lt.4095).and.
     &             (tof31(left,i,iadc).lt.3786)) THEN
         xhelp = tw31(left,i)/(tof31(left,i,iadc)**0.5)
         tof31(left,i,itdc) = tof31(left,i,itdc) + xhelp
         tdc_c(ch31a(i),hb31a(i))=tof31(left,i,itdc)
                                              ENDIF

         if ((tdc(ch31b(i),hb31b(i)).lt.4095).and.
     &             (tof31(right,i,iadc).lt.3786)) THEN
         xhelp = tw31(right,i)/(tof31(right,i,iadc)**0.5)
         tof31(right,i,itdc) = tof31(right,i,itdc) + xhelp
         tdc_c(ch31b(i),hb31b(i))=tof31(right,i,itdc)
                                             ENDIF
      ENDDO

      DO I=1,3
         if ((tdc(ch32a(i),hb32a(i)).lt.4095).and.
     &             (tof32(left,i,iadc).lt.3786)) THEN
         xhelp = tw32(left,i)/(tof32(left,i,iadc)**0.5)
         tof32(left,i,itdc) = tof32(left,i,itdc) + xhelp
         tdc_c(ch32a(i),hb32a(i))=tof32(left,i,itdc)
                                              ENDIF

         if ((tdc(ch32b(i),hb32b(i)).lt.4095).and.
     &             (tof32(right,i,iadc).lt.3786)) THEN
         xhelp = tw32(right,i)/(tof32(right,i,iadc)**0.5)
         tof32(right,i,itdc) = tof32(right,i,itdc) + xhelp
         tdc_c(ch32b(i),hb32b(i))=tof32(right,i,itdc)
                                             ENDIF
      ENDDO
      

C---------------------------------------------------------------
C--- calculate track position in paddle using timing difference
C--- now using the time-walk corrected TDC values
C---------------------------------------------------------------

      do i=1,3
         xtofpos(i)=100.
         ytofpos(i)=100.
      enddo 

C-----------------------------S1 --------------------------------

      IF (tof11_i.GT.none_find) THEN
         ytofpos(1)  = ((tof11(1,tof11_i,itdc)-tof11(2,tof11_i,itdc))/2.
     +        -y_coor_lin11(tof11_i,offset))/y_coor_lin11(tof11_i,slope)
       i=tof11_i
      endif

      IF (tof12_i.GT.none_find) THEN
         xtofpos(1)  = ((tof12(1,tof12_i,itdc)-tof12(2,tof12_i,itdc))/2.
     +        -x_coor_lin12(tof12_i,offset))/x_coor_lin12(tof12_i,slope)
      i=tof12_i
      endif


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

      IF (tof21_i.GT.none_find) THEN
         xtofpos(2) = ((tof21(1,tof21_i,itdc)-tof21(2,tof21_i,itdc))/2.
     +        -x_coor_lin21(tof21_i,offset))/x_coor_lin21(tof21_i,slope)
      i=tof21_i
      endif

      IF (tof22_i.GT.none_find) THEN
         ytofpos(2) = ((tof22(1,tof22_i,itdc)-tof22(2,tof22_i,itdc))/2.
     +        -y_coor_lin22(tof22_i,offset))/y_coor_lin22(tof22_i,slope)
      i=tof22_i
      endif
      

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

      IF (tof31_i.GT.none_find) THEN
         ytofpos(3)  = ((tof31(1,tof31_i,itdc)-tof31(2,tof31_i,itdc))/2.
     +        -y_coor_lin31(tof31_i,offset))/y_coor_lin31(tof31_i,slope)
      i=tof31_i
      endif

      IF (tof32_i.GT.none_find) THEN
         xtofpos(3)  = ((tof32(1,tof32_i,itdc)-tof32(2,tof32_i,itdc))/2.
     +        -x_coor_lin32(tof32_i,offset))/x_coor_lin32(tof32_i,slope)
      i=tof32_i
      endif


c      do i=1,3
c         if (abs(xtofpos(i)).gt.100.) then
c            xtofpos(i)=101.
c         endif
c         if (abs(ytofpos(i)).gt.100.) then
c            ytofpos(i)=101.
c         endif
c      enddo 


C--  restrict TDC measurements to physical paddle dimensions +/- 10 cm
C--  this cut is now stronger than in the old versions

        if (abs(xtofpos(1)).gt.31.)  xtofpos(1)=101.
        if (abs(xtofpos(2)).gt.19.)  xtofpos(2)=101.
        if (abs(xtofpos(3)).gt.19.)  xtofpos(3)=101.

        if (abs(ytofpos(1)).gt.26.)  ytofpos(1)=101.
        if (abs(ytofpos(2)).gt.18.)  ytofpos(2)=101.
        if (abs(ytofpos(3)).gt.18.)  ytofpos(3)=101.

C----------------------------------------------------------------------
C---------------------  zenith angle theta  ---------------------------
C----------------------------------------------------------------------
C-------------------  improved calculation  ---------------------------

       xhelp1=0.
       if (tof11_i.GT.none_find) xhelp1=tof11_x(tof11_i)
       if (xtofpos(1).lt.100)  xhelp1=xtofpos(1)

       yhelp1=0.
       if (tof12_i.GT.none_find) yhelp1=tof12_y(tof12_i)
       if (ytofpos(1).lt.100)  yhelp1=ytofpos(1)

       yhelp2=0.
       if (tof32_i.GT.none_find) yhelp2=tof32_y(tof32_i)
       if (ytofpos(3).lt.100)  yhelp2=ytofpos(3)

       xhelp2=0.
       if (tof31_i.GT.none_find) xhelp2=tof31_x(tof31_i)
       if (xtofpos(3).lt.100)  xhelp2=xtofpos(3)


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

       dx  = xhelp1 - xhelp2
       dy  = yhelp1 - yhelp2
       dr = sqrt(dx*dx+dy*dy)
       theta13 = atan(dr/tofarm13)


C------------------------------------------------------------------
C------------------------------------------------------------------
C-------angle and ADC(x) correction: moved to new dEdx routine
C------------------------------------------------------------------
C------------------------------------------------------------------

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

      dist = ZTOF(1) - ZTOF(5)
      c2 = (2.*0.01*dist)/(3.E08*50.E-12  )

       IF ((tof11_i.GT.none_find).AND.(tof31_i.GT.none_find).AND.
     &    (ytofpos(1).NE.101.).AND.(ytofpos(3).NE.101.)) 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
         if (iz.le.1) c1 = k_S11S31(1,ihelp)
         if (iz.eq.2) c1 = k_S11S31(2,ihelp)
         if (iz.gt.2) c1 = k_S11S31(3,ihelp)
         betatof_a(1) = c2/(cos(theta13)*(ds-c1))

C------- ToF Mask - S11 - S31

         tofmask(ch11a(tof11_i),hb11a(tof11_i)) =
     $        tofmask(ch11a(tof11_i),hb11a(tof11_i)) + 1
         tofmask(ch11b(tof11_i),hb11b(tof11_i)) =
     $        tofmask(ch11b(tof11_i),hb11b(tof11_i)) + 1

         tofmask(ch31a(tof31_i),hb31a(tof31_i)) =
     $        tofmask(ch31a(tof31_i),hb31a(tof31_i)) + 1
         tofmask(ch31b(tof31_i),hb31b(tof31_i)) = 
     $        tofmask(ch31b(tof31_i),hb31b(tof31_i)) + 1

C------- 
      
      ENDIF
      
C     S11 - S32

      dist = ZTOF(1) - ZTOF(6)
      c2 = (2.*0.01*dist)/(3.E08*50.E-12  )

       IF ((tof11_i.GT.none_find).AND.(tof32_i.GT.none_find).AND.
     &    (ytofpos(1).NE.101.).AND.(xtofpos(3).NE.101.)) 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
         if (iz.le.1) c1 = k_S11S32(1,ihelp)
         if (iz.eq.2) c1 = k_S11S32(2,ihelp)
         if (iz.gt.2) c1 = k_S11S32(3,ihelp)
         betatof_a(2) = c2/(cos(theta13)*(ds-c1))

C------- ToF Mask - S11 - S32

         tofmask(ch11a(tof11_i),hb11a(tof11_i)) =
     $        tofmask(ch11a(tof11_i),hb11a(tof11_i)) + 1
         tofmask(ch11b(tof11_i),hb11b(tof11_i)) =
     $        tofmask(ch11b(tof11_i),hb11b(tof11_i)) + 1

         tofmask(ch32a(tof32_i),hb32a(tof32_i)) =
     $        tofmask(ch32a(tof32_i),hb32a(tof32_i)) + 1
         tofmask(ch32b(tof32_i),hb32b(tof32_i)) = 
     $        tofmask(ch32b(tof32_i),hb32b(tof32_i)) + 1

C------- 

      ENDIF
      
C     S12 - S31

      dist = ZTOF(2) - ZTOF(5)
      c2 = (2.*0.01*dist)/(3.E08*50.E-12  )

       IF ((tof12_i.GT.none_find).AND.(tof31_i.GT.none_find).AND.
     &    (xtofpos(1).NE.101.).AND.(ytofpos(3).NE.101.)) 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
         if (iz.le.1) c1 = k_S12S31(1,ihelp)
         if (iz.eq.2) c1 = k_S12S31(2,ihelp)
         if (iz.gt.2) c1 = k_S12S31(3,ihelp)
         betatof_a(3) = c2/(cos(theta13)*(ds-c1))

C------- ToF Mask - S12 - S31

         tofmask(ch12a(tof12_i),hb12a(tof12_i)) =
     $        tofmask(ch12a(tof12_i),hb12a(tof12_i)) + 1
         tofmask(ch12b(tof12_i),hb12b(tof12_i)) =
     $        tofmask(ch12b(tof12_i),hb12b(tof12_i)) + 1

         tofmask(ch31a(tof31_i),hb31a(tof31_i)) =
     $        tofmask(ch31a(tof31_i),hb31a(tof31_i)) + 1
         tofmask(ch31b(tof31_i),hb31b(tof31_i)) = 
     $        tofmask(ch31b(tof31_i),hb31b(tof31_i)) + 1

C------- 

      ENDIF
            
C     S12 - S32

      dist = ZTOF(2) - ZTOF(6)
      c2 = (2.*0.01*dist)/(3.E08*50.E-12  )

       IF ((tof12_i.GT.none_find).AND.(tof32_i.GT.none_find).AND.
     &    (xtofpos(1).NE.101.).AND.(xtofpos(3).NE.101.)) 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
         if (iz.le.1) c1 = k_S12S32(1,ihelp)
         if (iz.eq.2) c1 = k_S12S32(2,ihelp)
         if (iz.gt.2) c1 = k_S12S32(3,ihelp)
         betatof_a(4) = c2/(cos(theta13)*(ds-c1))

C------- ToF Mask - S12 - S32

         tofmask(ch12a(tof12_i),hb12a(tof12_i)) =
     $        tofmask(ch12a(tof12_i),hb12a(tof12_i)) + 1
         tofmask(ch12b(tof12_i),hb12b(tof12_i)) =
     $        tofmask(ch12b(tof12_i),hb12b(tof12_i)) + 1

         tofmask(ch32a(tof32_i),hb32a(tof32_i)) =
     $        tofmask(ch32a(tof32_i),hb32a(tof32_i)) + 1
         tofmask(ch32b(tof32_i),hb32b(tof32_i)) = 
     $        tofmask(ch32b(tof32_i),hb32b(tof32_i)) + 1

C------- 

      ENDIF

C     S21 - S31

      dist = ZTOF(3) - ZTOF(5)
      c2 = (2.*0.01*dist)/(3.E08*50.E-12  )

       IF ((tof21_i.GT.none_find).AND.(tof31_i.GT.none_find).AND.
     &    (xtofpos(2).NE.101.).AND.(ytofpos(3).NE.101.)) 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
         if (iz.le.1) c1 = k_S21S31(1,ihelp)
         if (iz.eq.2) c1 = k_S21S31(2,ihelp)
         if (iz.gt.2) c1 = k_S21S31(3,ihelp)
         betatof_a(5) = c2/(cos(theta13)*(ds-c1))

C------- ToF Mask - S21 - S31

         tofmask(ch21a(tof21_i),hb21a(tof21_i)) =
     $        tofmask(ch21a(tof21_i),hb21a(tof21_i)) + 1
         tofmask(ch21b(tof21_i),hb21b(tof21_i)) =
     $        tofmask(ch21b(tof21_i),hb21b(tof21_i)) + 1

         tofmask(ch31a(tof31_i),hb31a(tof31_i)) =
     $        tofmask(ch31a(tof31_i),hb31a(tof31_i)) + 1
         tofmask(ch31b(tof31_i),hb31b(tof31_i)) = 
     $        tofmask(ch31b(tof31_i),hb31b(tof31_i)) + 1

C------- 

      ENDIF

C     S21 - S32

      dist = ZTOF(3) - ZTOF(6)
      c2 = (2.*0.01*dist)/(3.E08*50.E-12  )


       IF ((tof21_i.GT.none_find).AND.(tof32_i.GT.none_find).AND.
     &    (xtofpos(2).NE.101.).AND.(xtofpos(3).NE.101.)) 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
         if (iz.le.1) c1 = k_S21S32(1,ihelp)
         if (iz.eq.2) c1 = k_S21S32(2,ihelp)
         if (iz.gt.2) c1 = k_S21S32(3,ihelp)
         betatof_a(6) = c2/(cos(theta13)*(ds-c1))
                  
C------- ToF Mask - S21 - S32

         tofmask(ch21a(tof21_i),hb21a(tof21_i)) =
     $        tofmask(ch21a(tof21_i),hb21a(tof21_i)) + 1
         tofmask(ch21b(tof21_i),hb21b(tof21_i)) =
     $        tofmask(ch21b(tof21_i),hb21b(tof21_i)) + 1

         tofmask(ch32a(tof32_i),hb32a(tof32_i)) =
     $        tofmask(ch32a(tof32_i),hb32a(tof32_i)) + 1
         tofmask(ch32b(tof32_i),hb32b(tof32_i)) = 
     $        tofmask(ch32b(tof32_i),hb32b(tof32_i)) + 1

C------- 

      ENDIF

C     S22 - S31

      dist = ZTOF(4) - ZTOF(5)
      c2 = (2.*0.01*dist)/(3.E08*50.E-12  )

       IF ((tof22_i.GT.none_find).AND.(tof31_i.GT.none_find).AND.
     &    (ytofpos(2).NE.101.).AND.(ytofpos(3).NE.101.)) 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
         if (iz.le.1) c1 = k_S22S31(1,ihelp)
         if (iz.eq.2) c1 = k_S22S31(2,ihelp)
         if (iz.gt.2) c1 = k_S22S31(3,ihelp)
         betatof_a(7) = c2/(cos(theta13)*(ds-c1))

C------- ToF Mask - S22 - S31

         tofmask(ch22a(tof22_i),hb22a(tof22_i)) =
     $        tofmask(ch22a(tof22_i),hb22a(tof22_i)) + 1
         tofmask(ch22b(tof22_i),hb22b(tof22_i)) =
     $        tofmask(ch22b(tof22_i),hb22b(tof22_i)) + 1

         tofmask(ch31a(tof31_i),hb31a(tof31_i)) =
     $        tofmask(ch31a(tof31_i),hb31a(tof31_i)) + 1
         tofmask(ch31b(tof31_i),hb31b(tof31_i)) = 
     $        tofmask(ch31b(tof31_i),hb31b(tof31_i)) + 1

C-------          

      ENDIF
      
C     S22 - S32

      dist = ZTOF(4) - ZTOF(6)
      c2 = (2.*0.01*dist)/(3.E08*50.E-12  )

       IF ((tof22_i.GT.none_find).AND.(tof32_i.GT.none_find).AND.
     &    (ytofpos(2).NE.101.).AND.(xtofpos(3).NE.101.)) 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
         if (iz.le.1) c1 = k_S22S32(1,ihelp)
         if (iz.eq.2) c1 = k_S22S32(2,ihelp)
         if (iz.gt.2) c1 = k_S22S32(3,ihelp)
         betatof_a(8) = c2/(cos(theta13)*(ds-c1))

C------- ToF Mask - S22 - S32

         tofmask(ch22a(tof22_i),hb22a(tof22_i)) =
     $        tofmask(ch22a(tof22_i),hb22a(tof22_i)) + 1
         tofmask(ch22b(tof22_i),hb22b(tof22_i)) =
     $        tofmask(ch22b(tof22_i),hb22b(tof22_i)) + 1

         tofmask(ch32a(tof32_i),hb32a(tof32_i)) =
     $        tofmask(ch32a(tof32_i),hb32a(tof32_i)) + 1
         tofmask(ch32b(tof32_i),hb32b(tof32_i)) = 
     $        tofmask(ch32b(tof32_i),hb32b(tof32_i)) + 1

C-------   

      ENDIF

C     S11 - S21

      dist = ZTOF(1) - ZTOF(3)
      c2 = (2.*0.01*dist)/(3.E08*50.E-12  )

       IF ((tof11_i.GT.none_find).AND.(tof21_i.GT.none_find).AND.
     &    (ytofpos(1).NE.101.).AND.(xtofpos(2).NE.101.)) 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
         if (iz.le.1) c1 = k_S11S21(1,ihelp)
         if (iz.eq.2) c1 = k_S11S21(2,ihelp)
         if (iz.gt.2) c1 = k_S11S21(3,ihelp)
         betatof_a(9) = c2/(cos(theta13)*(ds-c1))

C------- ToF Mask - S11 - S21

         tofmask(ch11a(tof11_i),hb11a(tof11_i)) =
     $        tofmask(ch11a(tof11_i),hb11a(tof11_i)) + 1
         tofmask(ch11b(tof11_i),hb11b(tof11_i)) =
     $        tofmask(ch11b(tof11_i),hb11b(tof11_i)) + 1

         tofmask(ch21a(tof21_i),hb21a(tof21_i)) =
     $        tofmask(ch21a(tof21_i),hb21a(tof21_i)) + 1
         tofmask(ch21b(tof21_i),hb21b(tof21_i)) = 
     $        tofmask(ch21b(tof21_i),hb21b(tof21_i)) + 1

C-------   

      ENDIF
      
C     S11 - S22

      dist = ZTOF(1) - ZTOF(4)
      c2 = (2.*0.01*dist)/(3.E08*50.E-12  )

       IF ((tof11_i.GT.none_find).AND.(tof22_i.GT.none_find).AND.
     &    (ytofpos(1).NE.101.).AND.(ytofpos(2).NE.101.)) 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
         if (iz.le.1) c1 = k_S11S22(1,ihelp)
         if (iz.eq.2) c1 = k_S11S22(2,ihelp)
         if (iz.gt.2) c1 = k_S11S22(3,ihelp)
         betatof_a(10) = c2/(cos(theta13)*(ds-c1))

C------- ToF Mask - S11 - S22

         tofmask(ch11a(tof11_i),hb11a(tof11_i)) =
     $        tofmask(ch11a(tof11_i),hb11a(tof11_i)) + 1
         tofmask(ch11b(tof11_i),hb11b(tof11_i)) =
     $        tofmask(ch11b(tof11_i),hb11b(tof11_i)) + 1

         tofmask(ch22a(tof22_i),hb22a(tof22_i)) =
     $        tofmask(ch22a(tof22_i),hb22a(tof22_i)) + 1
         tofmask(ch22b(tof22_i),hb22b(tof22_i)) = 
     $        tofmask(ch22b(tof22_i),hb22b(tof22_i)) + 1

C-------   

      ENDIF

C     S12 - S21

      dist = ZTOF(2) - ZTOF(3)
      c2 = (2.*0.01*dist)/(3.E08*50.E-12  )

       IF ((tof12_i.GT.none_find).AND.(tof21_i.GT.none_find).AND.
     &    (xtofpos(1).NE.101.).AND.(xtofpos(2).NE.101.)) 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
         if (iz.le.1) c1 = k_S12S21(1,ihelp)
         if (iz.eq.2) c1 = k_S12S21(2,ihelp)
         if (iz.gt.2) c1 = k_S12S21(3,ihelp)
         betatof_a(11) = c2/(cos(theta13)*(ds-c1))

C------- ToF Mask - S12 - S21

         tofmask(ch12a(tof12_i),hb12a(tof12_i)) =
     $        tofmask(ch12a(tof12_i),hb12a(tof12_i)) + 1
         tofmask(ch12b(tof12_i),hb12b(tof12_i)) =
     $        tofmask(ch12b(tof12_i),hb12b(tof12_i)) + 1

         tofmask(ch21a(tof21_i),hb21a(tof21_i)) =
     $        tofmask(ch21a(tof21_i),hb21a(tof21_i)) + 1
         tofmask(ch21b(tof21_i),hb21b(tof21_i)) = 
     $        tofmask(ch21b(tof21_i),hb21b(tof21_i)) + 1

C-------  

      ENDIF

C     S12 - S22

      dist = ZTOF(2) - ZTOF(4)
      c2 = (2.*0.01*dist)/(3.E08*50.E-12  )

       IF ((tof12_i.GT.none_find).AND.(tof22_i.GT.none_find).AND.
     &    (xtofpos(1).NE.101.).AND.(ytofpos(2).NE.101.)) 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
         if (iz.le.1) c1 = k_S12S22(1,ihelp)
         if (iz.eq.2) c1 = k_S12S22(2,ihelp)
         if (iz.gt.2) c1 = k_S12S22(3,ihelp)
         betatof_a(12) = c2/(cos(theta13)*(ds-c1))

C------- ToF Mask - S12 - S22

         tofmask(ch12a(tof12_i),hb12a(tof12_i)) =
     $        tofmask(ch12a(tof12_i),hb12a(tof12_i)) + 1
         tofmask(ch12b(tof12_i),hb12b(tof12_i)) =
     $        tofmask(ch12b(tof12_i),hb12b(tof12_i)) + 1

         tofmask(ch22a(tof22_i),hb22a(tof22_i)) =
     $        tofmask(ch22a(tof22_i),hb22a(tof22_i)) + 1
         tofmask(ch22b(tof22_i),hb22b(tof22_i)) = 
     $        tofmask(ch22b(tof22_i),hb22b(tof22_i)) + 1

C-------  

      ENDIF

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

C--------  New mean beta  calculation  -----------------------

        do i=1,12
         btemp(i) = betatof_a(i)
        enddo

        betatof_a(13)=newbeta(1,btemp,hitvec,10.,10.,20.)

C--------------------------------------------------------------
C      write(*,*) betatof_a
c      write(*,*) xtofpos
c      write(*,*) ytofpos
c      write(*,*)'tofl2com beta', betatof_a
C      write(*,*) adcflagtof
c      write(*,*) 'tofl2com'
c      write(*,*) xtofpos
c      write(*,*) ytofpos
c      write(*,*) xtr_tof
c      write(*,*) ytr_tof
      
c 100  continue
        continue

C
      RETURN
      END


C------------------------------------------------------------------
C------------------------------------------------------------------

       function atten(is,ilay,ipad,x) 
       include  'input_tof.txt'
       real atten
       real x
       real xmin,xmax
       integer  ilay,ipad

*  S11 8 paddles  33.0 x 5.1 cm
*  S12 6 paddles  40.8 x 5.5 cm
*  S21 2 paddles  18.0 x 7.5 cm
*  S22 2 paddles  15.0 x 9.0 cm
*  S31 3 paddles  15.0 x 6.0 cm
*  S32 3 paddles  18.0 x 5.0 cm


c       if (ilay.eq.11) write(*,*) 'start ',ipad,is,adcx11(is,ipad,1),
c     &  adcx11(is,ipad,2),adcx11(is,ipad,3),adcx11(is,ipad,4)
c       if (ilay.eq.12) write(*,*) 'start ',ipad,is,adcx12(is,ipad,1),
c     &  adcx12(is,ipad,2),adcx12(is,ipad,3),adcx12(is,ipad,4)


       if (ilay.eq.11)  xmin=-33.0/2. 
       if (ilay.eq.11)  xmax= 33.0/2. 
       if (ilay.eq.12)  xmin=-40.8/2. 
       if (ilay.eq.12)  xmax= 40.8/2. 

       if (ilay.eq.21)  xmin=-18.0/2. 
       if (ilay.eq.21)  xmax= 18.0/2. 
       if (ilay.eq.22)  xmin=-15.0/2. 
       if (ilay.eq.22)  xmax= 15.0/2. 

       if (ilay.eq.31)  xmin=-15.0/2. 
       if (ilay.eq.31)  xmax= 15.0/2. 
       if (ilay.eq.32)  xmin=-18.0/2. 
       if (ilay.eq.32)  xmax= 18.0/2. 

       if (x .lt. xmin) x=xmin
       if (x .gt. xmax) x=xmax
 

       if (ilay.eq.11) atten=
     &    adcx11(is,ipad,1)*exp(x*adcx11(is,ipad,2))
     &  + adcx11(is,ipad,3)*exp(x*adcx11(is,ipad,4))

       if (ilay.eq.12) atten=
     &    adcx12(is,ipad,1)*exp(x*adcx12(is,ipad,2))
     &  + adcx12(is,ipad,3)*exp(x*adcx12(is,ipad,4))

       if (ilay.eq.21) atten=
     &    adcx21(is,ipad,1)*exp(x*adcx21(is,ipad,2))
     &  + adcx21(is,ipad,3)*exp(x*adcx21(is,ipad,4))

       if (ilay.eq.22) atten=
     &    adcx22(is,ipad,1)*exp(x*adcx22(is,ipad,2))
     &  + adcx22(is,ipad,3)*exp(x*adcx22(is,ipad,4))

       if (ilay.eq.31) atten=
     &    adcx31(is,ipad,1)*exp(x*adcx31(is,ipad,2))
     &  + adcx31(is,ipad,3)*exp(x*adcx31(is,ipad,4))

       if (ilay.eq.32) atten=
     &    adcx32(is,ipad,1)*exp(x*adcx32(is,ipad,2))
     &  + adcx32(is,ipad,3)*exp(x*adcx32(is,ipad,4))

        if (atten.gt.10000) atten=10000.

       end

C------------------------------------------------------------------
C------------------------------------------------------------------

       function pc_adc(ix) 
       include  'input_tof.txt'
       real pc_adc
       integer ix

       pc_adc=28.0407 + 0.628929*ix
     &   - 5.80901e-05*ix*ix + 3.14092e-08*ix*ix*ix
c       write(*,*) ix,pc_adc
       end

C------------------------------------------------------------------
C------------------------------------------------------------------

        function check_charge(theta,hitvec)

        include  'input_tof.txt'
        include  'tofcomm.txt'

        real check_charge  
        integer hitvec(6)  
        REAL CHARGE, theta

C  upper and lower limits  for the helium selection
        REAL A_l(24),A_h(24)
        DATA A_l /200,190,300,210,220,200,210,60,60,120,220,
     &  120,160,50,300,200,120,250,350,300,350,250,280,300/
        DATA A_h /550,490,800,600,650,600,600,260,200,380,
     &  620,380,550,200,850,560,400,750,900,800,880,800,750,800/

C The k1 constants for the beta calculation, only for S1-S3
C k2 constant is taken to be the standard 2D/c
        REAL k1(84)
        DATA k1 /50,59.3296,28.4328,-26.0818,5.91253,-19.588,
     &   -9.26316,24.7544,2.32465,-50.5058,-15.3195,-39.1443,
     &   -91.2546,-58.6243,-84.5641,-63.1516,-32.2091,-58.3358,
     &   13.8084,45.5322,33.2416,-11.5313,51.3271,75,-14.1141,
     &   42.8466,15.1794,-63.6672,-6.07739,-32.164,-41.771,10.5274,
     &   -9.46096,-81.7404,-28.783,-52.7167,-127.394,-69.6166,
     &   -93.4655,-98.9543,-42.863,-67.8244,-19.3238,31.1221,8.7319,
     &   -43.1627,5.55573,-14.4078,-83.4466,-47.4647,-77.8379,
     &   -108.222,-75.986,-101.297,-96.0205,-63.1881,-90.1372,
     &   -22.7347,8.31409,-19.6912,-7.49008,23.6979,-1.66677,
     &   1.81556,34.4668,6.23693,-100,-59.5861,-90.9159,-141.639,
     &   -89.2521,-112.881,-130.199,-77.0357,-98.4632,-60.2086,
     &   -4.82097,-29.3705,-43.6469,10.5884,-9.31304,-35.3329,
     &   25.2514,25.6/


        REAL zin(6)
        DATA zin /53.74, 53.04, 23.94, 23.44, -23.49, -24.34/

        REAL  c1,c2,xhelp,xhelp1,xhelp2,ds,dist,F
        REAL  sw,sxw,beta_mean_tof,w_i
        INTEGER  ihelp
        INTEGER ipmt(4)
        REAL time(4),beta1(4)

        REAL  adca(48), tdca(48)

        REAL  a1,a2
        INTEGER jj

c get rid of warnings EMILIANO
        i = 0
        slope = 0
        offset = 0
        none_find = 0
        none_ev = 0
        adc_ev = 0
        tdc_ev = 0
        iadc = 0
        itdc = 0
        right = 0
        left = 0
        tof12_y(1) = tof12_y(1) 
        tof11_x(1) = tof11_x(1)
        tof21_y(1) = tof21_y(1)
        tof22_x(1) = tof22_x(1)
        tof32_y(1) = tof32_y(1) 
        tof31_x(1) = tof31_x(1)
c get rid of warnings

C-----------------------------------------------------------
C--- get data
C-----------------------------------------------------------

         do j=1,8
         ih = 1 + 0 +((j-1)*2)
         adca(ih)   = adc(ch11a(j),hb11a(j))
         adca(ih+1) = adc(ch11b(j),hb11b(j))
         tdca(ih)   = tdc(ch11a(j),hb11a(j))
         tdca(ih+1) = tdc(ch11b(j),hb11b(j))
         enddo

         do j=1,6
         ih = 1 + 16+((j-1)*2)
         adca(ih)   = adc(ch12a(j),hb12a(j))
         adca(ih+1) = adc(ch12b(j),hb12b(j))
         tdca(ih)   = tdc(ch12a(j),hb12a(j))
         tdca(ih+1) = tdc(ch12b(j),hb12b(j))
         enddo

         do j=1,2
         ih = 1 + 28+((j-1)*2)
         adca(ih)   = adc(ch21a(j),hb21a(j))
         adca(ih+1) = adc(ch21b(j),hb21b(j))
         tdca(ih)   = tdc(ch21a(j),hb21a(j))
         tdca(ih+1) = tdc(ch21b(j),hb21b(j))
         enddo

         do j=1,2
         ih = 1 + 32+((j-1)*2)
         adca(ih)   = adc(ch22a(j),hb22a(j))
         adca(ih+1) = adc(ch22b(j),hb22b(j))
         tdca(ih)   = tdc(ch22a(j),hb22a(j))
         tdca(ih+1) = tdc(ch22b(j),hb22b(j))
         enddo

         do j=1,3
         ih = 1 + 36+((j-1)*2)
         adca(ih)   = adc(ch31a(j),hb31a(j))
         adca(ih+1) = adc(ch31b(j),hb31b(j))
         tdca(ih)   = tdc(ch31a(j),hb31a(j))
         tdca(ih+1) = tdc(ch31b(j),hb31b(j))
         enddo

         do j=1,3
         ih = 1 + 42+((j-1)*2)
         adca(ih)   = adc(ch32a(j),hb32a(j))
         adca(ih+1) = adc(ch32b(j),hb32b(j))
         tdca(ih)   = tdc(ch32a(j),hb32a(j))
         tdca(ih+1) = tdc(ch32b(j),hb32b(j))
         enddo


c         write(*,*) adca
c         write(*,*) tdca


C============   calculate beta and select charge > Z=1   ===============

        ICHARGE=1

C find hitted paddle by looking for ADC values on both sides
C since we looking for Z>1 this gives decent results

        tof11_i = hitvec(1)-1
        tof12_i = hitvec(2)-1
        tof21_i = hitvec(3)-1
        tof22_i = hitvec(4)-1
        tof31_i = hitvec(5)-1
        tof32_i = hitvec(6)-1

c        write(*,*) ' in charge check' 
c        write(*,*) theta,tof11_i,tof12_i,tof21_i,tof22_i,tof31_i,tof32_i

C----------------------------------------------------------------

        beta_help=100.
        beta_mean_tof=100.

        do jj=1,4
          beta1(jj) = 100.
        enddo

C----------------------------------------------------------------
C---------  S1 - S3 ---------------------------------------------
C----------------------------------------------------------------

C---------  S11 - S31 -------------------------------------------

        if ((tof11_i.gt.-1).and.(tof31_i.gt.-1)) then

        dist = zin(1) - zin(5)
        c2 = (2.*0.01*dist)/(3.E08*50.E-12)
        F = 1./cos(theta)

        ipmt(1)   = (tof11_i)*2+1
        ipmt(2)   = (tof11_i)*2+2
        ipmt(3)   = 36+(tof31_i)*2+1
        ipmt(4)   = 36+(tof31_i)*2+2

c        write(*,*) ipmt

        do jj=1,4
           time(jj) = tdca(ipmt(jj))
        enddo

c        write(*,*) time

        if ((time(1).lt.4095).and.(time(2).lt.4095).and.
     &     (time(3).lt.4095).and.(time(4).lt.4095)) then
         xhelp1 = time(1) + time(2)
         xhelp2 = time(3) + time(4)
         ds = xhelp1-xhelp2
         ihelp=0+(tof11_i)*3+tof31_i
         c1 = k1(ihelp+1)
         beta1(1) = c2*F/(ds-c1);
                                                 endif
c         write(*,*) beta1(1)
         endif  ! tof_....


C---------  S11 - S32 -------------------------------------------

        if ((tof11_i.gt.-1).and.(tof32_i.gt.-1)) then

        dist = zin(1) - zin(6)
        c2 = (2.*0.01*dist)/(3.E08*50.E-12)
        F = 1./cos(theta)

        ipmt(1)   = (tof11_i)*2+1
        ipmt(2)   = (tof11_i)*2+2
        ipmt(3)   = 42+(tof32_i)*2+1
        ipmt(4)   = 42+(tof32_i)*2+2

        do jj=1,4
           time(jj) = tdca(ipmt(jj))
        enddo

        if ((time(1).lt.4095).and.(time(2).lt.4095).and.
     &     (time(3).lt.4095).and.(time(4).lt.4095)) then
         xhelp1 = time(1) + time(2)
         xhelp2 = time(3) + time(4)
         ds = xhelp1-xhelp2
         ihelp=24+(tof11_i)*3+tof32_i
         c1 = k1(ihelp+1)
         beta1(2) = c2*F/(ds-c1);
                                                 endif
         endif  ! tof_....


C---------  S12 - S31 -------------------------------------------

        if ((tof12_i.gt.-1).and.(tof31_i.gt.-1)) then

        dist = zin(2) - zin(5)
        c2 = (2.*0.01*dist)/(3.E08*50.E-12)
        F = 1./cos(theta)

        ipmt(1)   = 16+(tof12_i)*2+1
        ipmt(2)   = 16+(tof12_i)*2+2
        ipmt(3)   = 36+(tof31_i)*2+1
        ipmt(4)   = 36+(tof31_i)*2+2

        do jj=1,4
           time(jj) = tdca(ipmt(jj))
        enddo

        if ((time(1).lt.4095).and.(time(2).lt.4095).and.
     &     (time(3).lt.4095).and.(time(4).lt.4095)) then
         xhelp1 = time(1) + time(2)
         xhelp2 = time(3) + time(4)
         ds = xhelp1-xhelp2
         ihelp=48+(tof12_i)*3+tof31_i
         c1 = k1(ihelp+1)
         beta1(3) = c2*F/(ds-c1);
                                                 endif
         endif  ! tof_....


C---------  S12 - S32 -------------------------------------------

        if ((tof12_i.gt.-1).and.(tof32_i.gt.-1)) then

        dist = zin(2) - zin(6)
        c2 = (2.*0.01*dist)/(3.E08*50.E-12)
        F = 1./cos(theta)

        ipmt(1)   = 16+(tof12_i)*2+1
        ipmt(2)   = 16+(tof12_i)*2+2
        ipmt(3)   = 42+(tof32_i)*2+1
        ipmt(4)   = 42+(tof32_i)*2+2

        do jj=1,4
           time(jj) = tdca(ipmt(jj))
        enddo

        if ((time(1).lt.4095).and.(time(2).lt.4095).and.
     &     (time(3).lt.4095).and.(time(4).lt.4095)) then
         xhelp1 = time(1) + time(2)
         xhelp2 = time(3) + time(4)
         ds = xhelp1-xhelp2
         ihelp=56+(tof12_i)*3+tof32_i
         c1 = k1(ihelp+1)
         beta1(4) = c2*F/(ds-c1);
                                                 endif

         endif  ! tof_....

c         write(*,*) beta1
  
C----  calculate  beta mean, only downward going particles are interesting ----

         sw=0.
         sxw=0.
         beta_mean_tof=100.

        do jj=1,4
        if ((beta1(jj).gt.0.1).and.(beta1(jj).lt.2.0)) then
            w_i=1./(0.13*0.13)
            sxw=sxw + beta1(jj)*w_i
            sw =sw + w_i ;
                                                      endif
        enddo

        if (sw.gt.0) beta_mean_tof=sxw/sw;

C        write(*,*) 'beta_mean_tof ',beta_mean_tof

        beta_help = beta_mean_tof  !  pow(beta_mean_tof,1.0) gave best results

CCCCC        endif  !  if tof11_i > -1 && ...... beta calculation

C-----------------------  Select charge   --------------------------

       charge=0
 
       if ((beta_mean_tof.gt.0.2).and.(beta_mean_tof.lt.2.0)) then

        icount1=0
        icount2=0
        icount3=0

        do jj=0,23
        a1 = adca(2*jj+1)
        a2 = adca(2*jj+2)
        if ((a1.lt.4095).and.(a2.lt.4095)) then
        a1 = adca(2*jj+1)*cos(theta)
        a2 = adca(2*jj+2)*cos(theta)
        xhelp  = 100000.
        xhelp1 = 100000.
        xhelp = sqrt(a1*a2)  ! geometric mean
        xhelp1 = beta_help*xhelp
C if geometric mean multiplied by beta_help is inside/outside helium 
C limits, increase counter
       if (xhelp1.lt.A_l(jj+1))  icount1=icount1+1
       if ((xhelp1.gt.A_l(jj+1)).and.(xhelp1.lt.A_h(jj+1)))
     &                           icount2=icount2+1
       if (xhelp1.gt.A_h(jj+1))  icount3=icount3+1
                                            endif
        enddo


C  if more than three paddles see the same...

        if (icount1 .gt. 3) charge=1
        if (icount2 .gt. 3) charge=2
        if (icount3 .gt. 3) charge=3

                                                    endif  ! 0.2<beta<2.0
 
C  no beta found? Sum up geometric means of paddles and derive the mean...
 
       if (beta_mean_tof.eq.100.) then

       xhelp  = 0.
       icount = 0 

       if (tof11_i.gt.-1) then
       jj=tof11_i
       a1 = adca(0+2*jj+1)
       a2 = adca(0+2*jj+2)
       if ((a1.lt.4095).and.(a2.lt.4095)) then
       a1 = a1*cos(theta)
       a2 = a2*cos(theta)
       xhelp = xhelp + sqrt(a1*a2)
       icount=icount+1
                    endif
                    endif

       if (tof12_i.gt.-1) then
       jj=tof12_i
       a1 = adca(16+2*jj+1)
       a2 = adca(16+2*jj+2)
       if ((a1.lt.4095).and.(a2.lt.4095)) then
       a1 = a1*cos(theta)
       a2 = a2*cos(theta)
       xhelp = xhelp + sqrt(a1*a2)
       icount=icount+1
                    endif
                    endif

       if (tof21_i.gt.-1) then
       jj=tof21_i
       a1 = adca(28+2*jj+1)
       a2 = adca(28+2*jj+2)
       if ((a1.lt.4095).and.(a2.lt.4095)) then
       a1 = a1*cos(theta)
       a2 = a2*cos(theta)
       xhelp = xhelp + sqrt(a1*a2)
       icount=icount+1
                    endif
                    endif

       if (tof22_i.gt.-1) then
       jj=tof22_i
       a1 = adca(32+2*jj+1)
       a2 = adca(32+2*jj+2)
       if ((a1.lt.4095).and.(a2.lt.4095)) then
       a1 = a1*cos(theta)
       a2 = a2*cos(theta)
       xhelp = xhelp + sqrt(a1*a2)
       icount=icount+1
                    endif
                    endif

       if (tof31_i.gt.-1) then
       jj=tof31_i
       a1 = adca(36+2*jj+1)
       a2 = adca(36+2*jj+2)
       if ((a1.lt.4095).and.(a2.lt.4095)) then
       a1 = a1*cos(theta)
       a2 = a2*cos(theta)
       xhelp = xhelp + sqrt(a1*a2)
       icount=icount+1
                    endif
                    endif

       if (tof32_i.gt.-1) then
       jj=tof32_i
       a1 = adca(42+2*jj+1)
       a2 = adca(42+2*jj+2)
       if ((a1.lt.4095).and.(a2.lt.4095)) then
       a1 = a1*cos(theta)
       a2 = a2*cos(theta)
       xhelp = xhelp + sqrt(a1*a2)
       icount=icount+1
                    endif
                    endif


       if (icount.gt.0) xhelp=xhelp/icount 
       if ((icount.gt.2).and.(xhelp.gt.1500.)) charge=3

                                  endif  ! beta_mean_tof.eq.100.

C        write(*,*) 'in function charge: ',beta_mean_tof,charge

        check_charge = charge


        END

C****************************************************************************
C****************************************************************************
C****************************************************************************

        function newbeta(iflag,b,hitvec,resmax,qualitycut,chi2cut)

        include  'input_tof.txt'
        include  'output_tof.txt'
        include  'tofcomm.txt'

        REAL newbeta 
        REAL resmax,qualitycut,chi2cut
        REAL w_i(12),w_il(6),quality,res,betachi,beta_mean_inv
        REAL sw,sxw,b(12),beta_mean,chi2,xhelp
c        REAL tdcfl(4,12)
        INTEGER tdcfl(4,12) !EM GCC4.7

        INTEGER iflag,icount,hitvec(6)

        INTEGER itop(12),ibot(12)
        DATA itop /1,1,2,2,3,3,4,4,1,1,2,2/
        DATA ibot /5,6,5,6,5,6,5,6,3,4,3,4/


c get rid of warnings EMILIANO
        slope = 0
        offset = 0
        none_find = 0
        none_ev = 0
        adc_ev = 0
        tdc_ev = 0
        iadc = 0
        itdc = 0
        right = 0
        left = 0
        tof12_y(1) = tof12_y(1) 
        tof11_x(1) = tof11_x(1)
        tof21_y(1) = tof21_y(1)
        tof22_x(1) = tof22_x(1)
        tof32_y(1) = tof32_y(1) 
        tof31_x(1) = tof31_x(1)
c get rid of warnings

C====================================================================

        tof11_i = hitvec(1)
        tof12_i = hitvec(2)
        tof21_i = hitvec(3)
        tof22_i = hitvec(4)
        tof31_i = hitvec(5)
        tof32_i = hitvec(6)

         if (iflag.eq.1) then   ! call from tofl2com
         do i=1,4
         do j=1,12
          tdcfl(i,j) =  tdcflagtof(i,j) 
         enddo
         enddo
                        endif

         if (iflag.eq.2) then   ! call from toftrk
         do i=1,4
         do j=1,12
          tdcfl(i,j) =  tdcflag(i,j) 
         enddo
         enddo
                        endif


C---  Find out ToF layers with artificial TDC values    -------------

        do jj=1,6
        w_il(jj) = 1000.
        enddo

 
        if (tof11_i.gt.0) then
        if ((tofmask(ch11a(tof11_i),hb11a(tof11_i)).gt.0).or.
     &   (tofmask(ch11b(tof11_i),hb11b(tof11_i)).gt.0)) then
        w_il(1)=0
        i1=tdcfl(ch11a(tof11_i),hb11a(tof11_i)) 
        i2=tdcfl(ch11b(tof11_i),hb11b(tof11_i)) 
        if ((i1.eq.1).or.(i2.eq.1)) w_il(1) = 1  ! tdcflag
                                                      endif
                           endif

        if (tof12_i.gt.0) then
        if ((tofmask(ch12a(tof12_i),hb12a(tof12_i)).gt.0).or.
     &   (tofmask(ch12b(tof12_i),hb12b(tof12_i)).gt.0)) then
        w_il(2)=0
        i1=tdcfl(ch12a(tof12_i),hb12a(tof12_i)) 
        i2=tdcfl(ch12b(tof12_i),hb12b(tof12_i)) 
        if ((i1.eq.1).or.(i2.eq.1)) w_il(2) = 1  ! tdcflag
                                                      endif
                           endif

        if (tof21_i.gt.0) then
        if ((tofmask(ch21a(tof21_i),hb21a(tof21_i)).gt.0).or.
     &   (tofmask(ch21b(tof21_i),hb21b(tof21_i)).gt.0)) then
        w_il(3)=0
        i1=tdcfl(ch21a(tof21_i),hb21a(tof21_i)) 
        i2=tdcfl(ch21b(tof21_i),hb21b(tof21_i)) 
        if ((i1.eq.1).or.(i2.eq.1)) w_il(3) = 1  ! tdcflag
                                                      endif
                           endif

        if (tof22_i.gt.0) then
        if ((tofmask(ch22a(tof22_i),hb22a(tof22_i)).gt.0).or.
     &   (tofmask(ch22b(tof22_i),hb22b(tof22_i)).gt.0)) then
        w_il(4)=0
        i1=tdcfl(ch22a(tof22_i),hb22a(tof22_i)) 
        i2=tdcfl(ch22b(tof22_i),hb22b(tof22_i)) 
        if ((i1.eq.1).or.(i2.eq.1)) w_il(4) = 1  ! tdcflag
                                                      endif
                           endif

        if (tof31_i.gt.0) then
        if ((tofmask(ch31a(tof31_i),hb11a(tof31_i)).gt.0).or.
     &   (tofmask(ch31b(tof31_i),hb31b(tof31_i)).gt.0)) then
        w_il(5)=0
        i1=tdcfl(ch31a(tof31_i),hb31a(tof31_i)) 
        i2=tdcfl(ch31b(tof31_i),hb31b(tof31_i)) 
        if ((i1.eq.1).or.(i2.eq.1)) w_il(5) = 1  ! tdcflag
                                                      endif
                           endif

        if (tof32_i.gt.0) then
        if ((tofmask(ch32a(tof32_i),hb32a(tof32_i)).gt.0).or.
     &   (tofmask(ch32b(tof32_i),hb32b(tof32_i)).gt.0)) then
        w_il(6)=0
        i1=tdcfl(ch32a(tof32_i),hb32a(tof32_i)) 
        i2=tdcfl(ch32b(tof32_i),hb32b(tof32_i)) 
        if ((i1.eq.1).or.(i2.eq.1)) w_il(6) = 1  ! tdcflag
                                                      endif
                           endif

C------------------------------------------------------------------------
C---  Set weights for the 12 measurements using information for top and bottom:
C---  if no measurements: weight = set to very high value=> not used
C---  top or bottom artificial: weight*sqrt(2)
C---  top and bottom artificial: weight*sqrt(2)*sqrt(2)

       DO jj=1,12
       if (jj.le.4)           xhelp = 0.11            ! S1-S3
       if ((jj.gt.4).and.(jj.le.8)) xhelp = 0.18      ! S2-S3
       if (jj.gt.8)           xhelp = 0.28            ! S1-S2
       if ((w_il(itop(jj)).eq.1000.).and.(w_il(ibot(jj)).eq.1000.)) 
     &   xhelp = 1.E09
       if ((w_il(itop(jj)).eq.1).or.(w_il(ibot(jj)).eq.1.)) 
     &   xhelp = xhelp*1.414 
       if ((w_il(itop(jj)).eq.1).and.(w_il(ibot(jj)).eq.1.)) 
     &   xhelp = xhelp*2. 
       w_i(jj) = 1./xhelp
       ENDDO

C========================================================================
C--- Calculate mean beta for the first time -----------------------------
C--- We are using "1/beta" since its error is gaussian ------------------

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

        DO jj=1,12
        IF ((abs(1./b(jj)).gt.0.1).and.(abs(1./b(jj)).lt.15.)) THEN
            icount = icount+1
            sxw    = sxw + (1./b(jj))*w_i(jj)*w_i(jj) 
            sw     = sw + w_i(jj)*w_i(jj)
        ENDIF
        ENDDO
 
        if (icount.gt.0) beta_mean=1./(sxw/sw)
        beta_mean_inv = 1./beta_mean

  
C--- Calculate beta for the second time, use residuals of the single
C--- measurements to get a chi2 value

        icount  = 0
        sw      = 0.
        sxw     = 0.
        betachi = 100.
        chi2    = 0.
        quality = 0.

        DO jj=1,12
        IF ((abs(1./b(jj)).gt.0.1).and.(abs(1./b(jj)).lt.15.)
     &                                .and.(w_i(jj).GT.0.01)) THEN
            res    = beta_mean_inv - (1./b(jj)) ;
            if (abs(res*w_i(jj)).lt.resmax) THEN
            chi2   = chi2 + (res*w_i(jj))**2.
            icount = icount+1
            sxw    = sxw + (1./b(jj))*w_i(jj)*w_i(jj)
            sw     = sw + w_i(jj)*w_i(jj)
                                            ENDIF
        ENDIF
        ENDDO

c        quality = sw 
        quality = sqrt(sw) 

        if (icount.eq.0) chi2 = 1000.
        if (icount.gt.0) chi2 = chi2/(icount) 

        if (icount.gt.0) betachi=1./(sxw/sw);

        beta_mean=100.
        if ((chi2.lt.chi2cut).and.(quality.gt.qualitycut)) 
     &  beta_mean = betachi
        newbeta = beta_mean

        END

C****************************************************************************
