/[PAMELA software]/DarthVader/ToFLevel2/src/tofl2com.for

Diff of /DarthVader/ToFLevel2/src/tofl2com.for

Parent Directory | Revision Log | View Patch Patch

-revision 1.6 by mocchiut,
Mon Apr 30 15:46:30 2007 UTC
+revision 1.12 by mocchiut,
Tue Aug 11 12:58:24 2009 UTC
 Line 1
- ******************************************************************************
+ C******************************************************************************
- *
+ C
- *  08-12-06 WM: adc_c-bug :  The raw ADc value was multiplied with cos(theta)
+ C  08-12-06 WM: adc_c-bug :  The raw ADc value was multiplied with cos(theta)
- *  and AFTER that there was an if statement "if tof32(right,i,iadc) < 4095"
+ C  and AFTER that there was an if statement "if tof32(right,i,iadc) < 4095"
- *
+ C
- *  jan-07 GF: ADCflags(4,12) inserted to flag artificial ADC values
+ C  jan-07 GF: ADCflags(4,12) inserted to flag artificial ADC values
- *  jan-07 WM: artificial ADC values created using attenuation calibration
+ C  jan-07 WM: artificial ADC values created using attenuation calibration
- *  jan-07 WM: modified xtofpos flag "101". xtofpos must be inside physical
+ C  jan-07 WM: modified xtofpos flag "101". xtofpos must be inside physical
- *             dimension of the paddle +/- 10 cm
+ C             dimension of the paddle +/- 10 cm
- *  jan-07 WM: if xtofpos=101 then this paddle is not used for beta
+ C  jan-07 WM: if xtofpos=101 then this paddle is not used for beta
- *             calculation
+ C             calculation
- *  jan-07 WM: the definition for a "hit" is changed: Now we must have a
+ C  jan-07 WM: the definition for a "hit" is changed: Now we must have a
- *             valid TDC signal on both sides
+ C             valid TDC signal on both sides
- *  jan-07 WM: flag for PMTs #10 and #35 added, TDC=819 due to bit-shift
+ C  jan-07 WM: flag for PMTs #10 and #35 added, TDC=819 due to bit-shift
- *  jan-07 WM: bug removed: in some cases tdc_tw was calculated due to a
+ C  jan-07 WM: bug removed: in some cases tdc_tw was calculated due to a
- *             leftover "xhelp" value
+ C             leftover "xhelp" value
- *  apr-07 WM: attenuation fit curve is now a double exponential fit
+ C  apr-07 WM: attenuation fit curve is now a double exponential fit
- *             conversion from raw ADC to pC using calibration function
+ C             conversion from raw ADC to pC using calibration function
- *             variables xtr_tof and ytr_tof inserted (filled with default)
+ 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******************************************************************************
        INTEGER FUNCTION TOFL2COM()
  c
-Line 34 
 C
+Line 47 
 C
        LOGICAL check
        REAL secure
-       INTEGER j
+       INTEGER j,hitvec(6)
-       REAL xhelp_a,xhelp_t
        REAL dx,dy,dr,ds
-       REAL yhelp,xhelp,xhelp1,xhelp2
+       REAL yhelp,yhelp1,yhelp2,xhelp,xhelp1,xhelp2
-       REAL c1,c2,sw,sxw,w_i
+       REAL c1,c2
-       INTEGER icount
- c      REAL xdummy
+ 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
-       REAL beta_mean
  c     value for status of each PM-data
  c     first index  : 1 = left, 2 = right
  c     second index : 1... number of paddle
-Line 57 
 c     second index : 1... number of padd
+Line 67 
 c     second index : 1... number of padd
        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--   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
-Line 70 
 C--   DATA ZTOF/53.74,53.04,23.94,23.44,
+Line 115 
 C--   DATA ZTOF/53.74,53.04,23.94,23.44,
        REAL hepratio
        INTEGER ihelp
-       REAL xkorr
+       REAL xkorr,btemp(12)
+       REAL atten,pc_adc,check_charge,newbeta
+       INTEGER IZ
+       REAL k1corrA1,k1corrB1,k1corrC1
-       real atten,pc_adc
+       INTEGER ifst
+       DATA  ifst /0/
  C---------------------------------------
  C
-Line 82 
 C
+Line 134 
 C
  C
  C     CALCULATE COMMON VARIABLES
  C
+ C-------------------------------------------------------------------
- *******************************************************************
+         if (ifst.eq.0) then
-       icounter = icounter + 1
- *     amplitude has to be 'secure' higher than pedestal for an adc event
+         ifst=1
-       secure = 2.
+ 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.  !
+        hepratio = 4.  !
-       offset = 1
+        offset = 1
-       slope = 2
+        slope = 2
-       left = 1
+        left = 1
-       right = 2
+        right = 2
-       none_ev = 0
+        none_ev = 0
-       none_find = 0
+        none_find = 0
-       tdc_ev = 1
+        tdc_ev = 1
-       adc_ev = 1
+        adc_ev = 1
-       itdc = 1
+        itdc = 1
-       iadc = 2
+        iadc = 2
+ C--- These are the corrections to the k1-value for Z>2 particles
+        k1corrA1 = 0.
+        k1corrB1 = -5.0
+        k1corrC1=  8.0
+         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.
-Line 151 
 C---  since this is standalone ToF fill
+Line 221 
 C---  since this is standalone ToF fill
  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
-Line 483 
 c     check if an other paddle has also
+Line 552 
 c     check if an other paddle has also
                 ENDIF
              ENDIF
           ENDIF
-       ENDDO
+        ENDDO
-       do i=1,6
+        do i=1,6
           tof_i_flag(i)=0
           tof_j_flag(i)=0
-       enddo
+        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
-       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
  C------------------------------------------------------------------
- C---  calculate track position in paddle using timing difference
+ 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
+        do i=1,3
           xtofpos(i)=100.
           ytofpos(i)=100.
-       enddo
+        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)
+      +   -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_lin12(tof12_i,offset))/x_coor_lin12(tof12_i,slope)
+      +   -x_coor_lin12c(tof12_i,offset))/x_coor_lin12c(tof12_i,slope)
        endif
-Line 530 
 C-----------------------------S2 -------
+Line 611 
 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)
+      +    -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_lin22(tof22_i,offset))/y_coor_lin22(tof22_i,slope)
+      +    -y_coor_lin22c(tof22_i,offset))/y_coor_lin22c(tof22_i,slope)
        endif
-Line 543 
 C-----------------------------S3 -------
+Line 624 
 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)
+      +    -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_lin32(tof32_i,offset))/x_coor_lin32(tof32_i,slope)
+      +    -x_coor_lin32c(tof32_i,offset))/x_coor_lin32c(tof32_i,slope)
        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-----------------------------  old  ----------------------------------
-       dx=0.
-       dy=0.
-       dr=0.
-       theta13 = 0.
-          IF ((tof12_i.GT.none_find).AND.(tof32_i.GT.none_find))
-      &        dx  = xtofpos(1) - xtofpos(3)
-          IF ((tof11_i.GT.none_find).AND.(tof31_i.GT.none_find))
-      &        dy  = ytofpos(1) - ytofpos(3)
-          dr = sqrt(dx*dx+dy*dy)
-          theta13 = atan(dr/tofarm13)
- C------------------------------------------------------------------
  c      dx=0.
  c      dy=0.
  c      dr=0.
- c      theta12 = 0.
+ c      theta13 = 0.
  c
- c         IF ((tof12_i.GT.none_find).AND.(tof21_i.GT.none_find))
+ c         IF ((tof12_i.GT.none_find).AND.(tof32_i.GT.none_find))
- c     &        dx  = xtofpos(1) - xtofpos(2)
+ c     &        dx  = xtofpos(1) - xtofpos(3)
- c         IF ((tof11_i.GT.none_find).AND.(tof22_i.GT.none_find))
+ c         IF ((tof11_i.GT.none_find).AND.(tof31_i.GT.none_find))
- c     &        dy  = ytofpos(1) - ytofpos(2)
+ c     &        dy  = ytofpos(1) - ytofpos(3)
  c         dr = sqrt(dx*dx+dy*dy)
- c         theta12 = atan(dr/tofarm12)
+ c         theta13 = atan(dr/tofarm13)
- c
- c      dx=0.
- c      dy=0.
- c      dr=0.
- c      theta23 = 0.
  c
- c         IF ((tof21_i.GT.none_find).AND.(tof32_i.GT.none_find))
+ c
- c     &        dx  = xtofpos(2) - xtofpos(3)
+ C-----------------------------  new  ------------------------------
- c         IF ((tof22_i.GT.none_find).AND.(tof31_i.GT.none_find))
- c     &        dy  = ytofpos(2) - ytofpos(3)
- c         dr = sqrt(dx*dx+dy*dy)
- c         theta23 = atan(dr/tofarm23)
- c
- 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
-Line 632 
 c       DATA tof32_y/ -5.0,0.0,5.0/
+Line 707 
 c       DATA tof32_y/ -5.0,0.0,5.0/
  C----------------------------  S1 -------------------------------------
-        yhelp=0.
+ 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
- c         if (tof11(left,i,iadc).eq.4095) then
           if (adc(ch11a(i),hb11a(i)).eq.4095) then
              xkorr = atten(left,11,i,yhelp)
-             xkorr=xkorr/hepratio
+             if (iz.le.1) xkorr=xkorr/hepratio
              tof11(left,i,iadc)=xkorr/cos(theta13)
- c            write(*,*) 'tofl2 left ',i, tof11(left,i,iadc)
              adcflagtof(ch11a(i),hb11a(i)) = 1
           endif
- c         if (tof11(right,i,iadc).eq.4095) then
           if (adc(ch11b(i),hb11b(i)).eq.4095) then
              xkorr = atten(right,11,i,yhelp)
-             xkorr=xkorr/hepratio
+             if (iz.le.1) xkorr=xkorr/hepratio
              tof11(right,i,iadc)=xkorr/cos(theta13)
- c            write(*,*) 'tofl2 right ',i, tof11(right,i,iadc)
              adcflagtof(ch11b(i),hb11b(i)) = 1
           endif
         ENDIF
-        xhelp=0.
+ 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
- c         if (tof12(left,i,iadc).eq.4095) then
           if (adc(ch12a(i),hb12a(i)).eq.4095) then
              xkorr = atten(left,12,i,xhelp)
-             xkorr=xkorr/hepratio
+             if (iz.le.1) xkorr=xkorr/hepratio
              tof12(left,i,iadc) = xkorr/cos(theta13)
              adcflagtof(ch12a(i),hb12a(i)) = 1
           endif
- c         if (tof12(right,i,iadc).eq.4095) then
           if (adc(ch12b(i),hb12b(i)).eq.4095) then
              xkorr = atten(right,12,i,xhelp)
-             xkorr=xkorr/hepratio
+             if (iz.le.1) xkorr=xkorr/hepratio
              tof12(right,i,iadc) = xkorr/cos(theta13)
              adcflagtof(ch12b(i),hb12b(i)) = 1
           endif
-Line 680 
 c         if (tof12(right,i,iadc).eq.409
+Line 751 
 c         if (tof12(right,i,iadc).eq.409
  C-----------------------------S2 --------------------------------
-        xhelp=0.
+ 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
- c         if (tof21(left,i,iadc).eq.4095) then
           if (adc(ch21a(i),hb21a(i)).eq.4095) then
              xkorr = atten(left,21,i,xhelp)
-             xkorr=xkorr/hepratio
+             if (iz.le.1) xkorr=xkorr/hepratio
              tof21(left,i,iadc) = xkorr/cos(theta13)
              adcflagtof(ch21a(i),hb21a(i)) = 1
           endif
- c         if (tof21(right,i,iadc).eq.4095) then
           if (adc(ch21b(i),hb21b(i)).eq.4095) then
              xkorr = atten(right,21,i,xhelp)
-             xkorr=xkorr/hepratio
+             if (iz.le.1) xkorr=xkorr/hepratio
              tof21(right,i,iadc) = xkorr/cos(theta13)
              adcflagtof(ch21b(i),hb21b(i)) = 1
           endif
         ENDIF
-        yhelp=0.
+ 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
- c         if (tof22(left,i,iadc).eq.4095) then
           if (adc(ch22a(i),hb22a(i)).eq.4095) then
              xkorr = atten(left,22,i,yhelp)
-             xkorr=xkorr/hepratio
+             if (iz.le.1) xkorr=xkorr/hepratio
              tof22(left,i,iadc) = xkorr/cos(theta13)
              adcflagtof(ch22a(i),hb22a(i)) = 1
           endif
- c         if (tof22(right,i,iadc).eq.4095) then
           if (adc(ch22b(i),hb22b(i)).eq.4095) then
              xkorr = atten(right,22,i,yhelp)
-             xkorr=xkorr/hepratio
+             if (iz.le.1) xkorr=xkorr/hepratio
              tof22(right,i,iadc) = xkorr/cos(theta13)
              adcflagtof(ch22b(i),hb22b(i)) = 1
           endif
-Line 727 
 c         if (tof22(right,i,iadc).eq.409
+Line 796 
 c         if (tof22(right,i,iadc).eq.409
  C-----------------------------S3 --------------------------------
-        yhelp=0.
+ 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
- c         if (tof31(left,i,iadc).eq.4095) then
           if (adc(ch31a(i),hb31a(i)).eq.4095) then
              xkorr = atten(left,31,i,yhelp)
-             xkorr=xkorr/hepratio
+             if (iz.le.1) xkorr=xkorr/hepratio
              tof31(left,i,iadc) = xkorr/cos(theta13)
              adcflagtof(ch31a(i),hb31a(i)) = 1
           endif
- c         if (tof31(right,i,iadc).eq.4095) then
           if (adc(ch31b(i),hb31b(i)).eq.4095) then
              xkorr = atten(right,31,i,yhelp)
-             xkorr=xkorr/hepratio
+             if (iz.le.1) xkorr=xkorr/hepratio
              tof31(right,i,iadc) = xkorr/cos(theta13)
              adcflagtof(ch31b(i),hb31b(i)) = 1
           endif
         ENDIF
-        xhelp=0.
+ 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
- c         if (tof32(left,i,iadc).eq.4095) then
           if (adc(ch32a(i),hb32a(i)).eq.4095) then
              xkorr = atten(left,32,i,xhelp)
-             xkorr=xkorr/hepratio
+             if (iz.le.1) xkorr=xkorr/hepratio
              tof32(left,i,iadc) = xkorr/cos(theta13)
              adcflagtof(ch32a(i),hb32a(i)) = 1
           endif
- c         if (tof32(right,i,iadc).eq.4095) then
           if (adc(ch32b(i),hb32b(i)).eq.4095) then
              xkorr = atten(right,32,i,xhelp)
-             xkorr=xkorr/hepratio
+             if (iz.le.1) xkorr=xkorr/hepratio
              tof32(right,i,iadc) = xkorr/cos(theta13)
              adcflagtof(ch32b(i),hb32b(i)) = 1
           endif
         ENDIF
- C--------------------------------------------------------------------
+ C-------------------------------------------------------------------
  C--------------------Time walk correction  -------------------------
- C--------------------------------------------------------------------
+ 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
-        xhelp= 0.
+          if ((tdc(ch11a(i),hb11a(i)).lt.4095).and.
-        xhelp_a = tof11(left,i,iadc)
+      &             (tof11(left,i,iadc).lt.3786)) THEN
-        xhelp_t = tof11(left,i,itdc)
+          xhelp = tw11(left,i)/(tof11(left,i,iadc)**0.5)
- c       if (xhelp_a .eq.0) write (*,*) '11 ',i,xhelp_a
+          tof11(left,i,itdc) = tof11(left,i,itdc) + xhelp
-        if(xhelp_a<3786) xhelp = tw11(left,i)/sqrt(xhelp_a)
+          tdc_c(ch11a(i),hb11a(i))=tof11(left,i,itdc)
-        tof11(left,i,itdc) = xhelp_t  + xhelp
+                                               ENDIF
-        tdc_c(ch11a(i),hb11a(i))=tof11(left,i,itdc)
-        xhelp_a = tof11(right,i,iadc)
+          if ((tdc(ch11b(i),hb11b(i)).lt.4095).and.
-        xhelp_t = tof11(right,i,itdc)
+      &             (tof11(right,i,iadc).lt.3786)) THEN
-        if(xhelp_a<3786) xhelp = tw11(right,i)/sqrt(xhelp_a)
+          xhelp = tw11(right,i)/(tof11(right,i,iadc)**0.5)
-        tof11(right,i,itdc) = xhelp_t  + xhelp
+          tof11(right,i,itdc) = tof11(right,i,itdc) + xhelp
-        tdc_c(ch11b(i),hb11b(i))=tof11(right,i,itdc)
+          tdc_c(ch11b(i),hb11b(i))=tof11(right,i,itdc)
+                                              ENDIF
        ENDDO
        DO i=1,6
-        xhelp= 0.
+          if ((tdc(ch12a(i),hb12a(i)).lt.4095).and.
-        xhelp_a = tof12(left,i,iadc)
+      &             (tof12(left,i,iadc).lt.3786)) THEN
-        xhelp_t = tof12(left,i,itdc)
+          xhelp = tw12(left,i)/(tof12(left,i,iadc)**0.5)
- c       if (xhelp_a .eq.0) write (*,*) '12 ',i,xhelp_a
+          tof12(left,i,itdc) = tof12(left,i,itdc) + xhelp
-        if(xhelp_a<3786) xhelp = tw12(left,i)/sqrt(xhelp_a)
+          tdc_c(ch12a(i),hb12a(i))=tof12(left,i,itdc)
-        tof12(left,i,itdc) = xhelp_t  + xhelp
+                                               ENDIF
-        tdc_c(ch12a(i),hb12a(i))=tof12(left,i,itdc)
-        xhelp_a = tof12(right,i,iadc)
+          if ((tdc(ch12b(i),hb12b(i)).lt.4095).and.
-        xhelp_t = tof12(right,i,itdc)
+      &             (tof12(right,i,iadc).lt.3786)) THEN
-        if(xhelp_a<3786) xhelp = tw12(right,i)/sqrt(xhelp_a)
+          xhelp = tw12(right,i)/(tof12(right,i,iadc)**0.5)
-        tof12(right,i,itdc) = xhelp_t  + xhelp
+          tof12(right,i,itdc) = tof12(right,i,itdc) + xhelp
-        tdc_c(ch12b(i),hb12b(i))=tof12(right,i,itdc)
+          tdc_c(ch12b(i),hb12b(i))=tof12(right,i,itdc)
+                                              ENDIF
        ENDDO
  C----
-       DO i=1,2
+       DO I=1,2
-        xhelp= 0.
+          if ((tdc(ch21a(i),hb21a(i)).lt.4095).and.
-        xhelp_a = tof21(left,i,iadc)
+      &             (tof21(left,i,iadc).lt.3786)) THEN
-        xhelp_t = tof21(left,i,itdc)
+          xhelp = tw21(left,i)/(tof21(left,i,iadc)**0.5)
- c       if (xhelp_a .eq.0) write (*,*) '21 ',i,xhelp_a
+          tof21(left,i,itdc) = tof21(left,i,itdc) + xhelp
-        if(xhelp_a<3786) xhelp = tw21(left,i)/sqrt(xhelp_a)
+          tdc_c(ch21a(i),hb21a(i))=tof21(left,i,itdc)
-        tof21(left,i,itdc) = xhelp_t  + xhelp
+                                               ENDIF
-        tdc_c(ch21a(i),hb21a(i))=tof21(left,i,itdc)
-        xhelp_a = tof21(right,i,iadc)
+          if ((tdc(ch21b(i),hb21b(i)).lt.4095).and.
-        xhelp_t = tof21(right,i,itdc)
+      &             (tof21(right,i,iadc).lt.3786)) THEN
-        if(xhelp_a<3786) xhelp = tw21(right,i)/sqrt(xhelp_a)
+          xhelp = tw21(right,i)/(tof21(right,i,iadc)**0.5)
-        tof21(right,i,itdc) = xhelp_t  + xhelp
+          tof21(right,i,itdc) = tof21(right,i,itdc) + xhelp
-        tdc_c(ch21b(i),hb21b(i))=tof21(right,i,itdc)
+          tdc_c(ch21b(i),hb21b(i))=tof21(right,i,itdc)
-       ENDDO
+                                              ENDIF
+       ENDDO
-       DO i=1,2
-        xhelp= 0.
+       DO I=1,2
-        xhelp_a = tof22(left,i,iadc)
+          if ((tdc(ch22a(i),hb22a(i)).lt.4095).and.
-        xhelp_t = tof22(left,i,itdc)
+      &             (tof22(left,i,iadc).lt.3786)) THEN
- c       if (xhelp_a .eq.0) write (*,*) '22 ',i,xhelp_a
+          xhelp = tw22(left,i)/(tof22(left,i,iadc)**0.5)
+          tof22(left,i,itdc) = tof22(left,i,itdc) + xhelp
-        if(xhelp_a<3786) xhelp = tw22(left,i)/sqrt(xhelp_a)
+          tdc_c(ch22a(i),hb22a(i))=tof22(left,i,itdc)
-        tof22(left,i,itdc) = xhelp_t  + xhelp
+                                               ENDIF
-        tdc_c(ch22a(i),hb22a(i))=tof22(left,i,itdc)
-        xhelp_a = tof22(right,i,iadc)
+          if ((tdc(ch22b(i),hb22b(i)).lt.4095).and.
-        xhelp_t = tof22(right,i,itdc)
+      &             (tof22(right,i,iadc).lt.3786)) THEN
-        if(xhelp_a<3786) xhelp = tw22(right,i)/sqrt(xhelp_a)
+          xhelp = tw22(right,i)/(tof22(right,i,iadc)**0.5)
-        tof22(right,i,itdc) = xhelp_t  + xhelp
+          tof22(right,i,itdc) = tof22(right,i,itdc) + xhelp
-        tdc_c(ch22b(i),hb22b(i))=tof22(right,i,itdc)
+          tdc_c(ch22b(i),hb22b(i))=tof22(right,i,itdc)
+                                              ENDIF
        ENDDO
- C----
-       DO i=1,3
+ C----
-        xhelp= 0.
+       DO I=1,3
-        xhelp_a = tof31(left,i,iadc)
+          if ((tdc(ch31a(i),hb31a(i)).lt.4095).and.
-        xhelp_t = tof31(left,i,itdc)
+      &             (tof31(left,i,iadc).lt.3786)) THEN
- c       if (xhelp_a .eq.0) write (*,*) '31 ',i,xhelp_a
+          xhelp = tw31(left,i)/(tof31(left,i,iadc)**0.5)
+          tof31(left,i,itdc) = tof31(left,i,itdc) + xhelp
-        if(xhelp_a<3786) xhelp = tw31(left,i)/sqrt(xhelp_a)
+          tdc_c(ch31a(i),hb31a(i))=tof31(left,i,itdc)
-        tof31(left,i,itdc) = xhelp_t  + xhelp
+                                               ENDIF
-        tdc_c(ch31a(i),hb31a(i))=tof31(left,i,itdc)
-        xhelp_a = tof31(right,i,iadc)
+          if ((tdc(ch31b(i),hb31b(i)).lt.4095).and.
-        xhelp_t = tof31(right,i,itdc)
+      &             (tof31(right,i,iadc).lt.3786)) THEN
-        if(xhelp_a<3786) xhelp = tw31(right,i)/sqrt(xhelp_a)
+          xhelp = tw31(right,i)/(tof31(right,i,iadc)**0.5)
-        tof31(right,i,itdc) = xhelp_t  + xhelp
+          tof31(right,i,itdc) = tof31(right,i,itdc) + xhelp
-        tdc_c(ch31b(i),hb31b(i))=tof31(right,i,itdc)
+          tdc_c(ch31b(i),hb31b(i))=tof31(right,i,itdc)
-       ENDDO
+                                              ENDIF
+       ENDDO
-       DO i=1,3
-        xhelp= 0.
+       DO I=1,3
-        xhelp_a = tof32(left,i,iadc)
+          if ((tdc(ch32a(i),hb32a(i)).lt.4095).and.
-        xhelp_t = tof32(left,i,itdc)
+      &             (tof32(left,i,iadc).lt.3786)) THEN
- c       if (xhelp_a .eq.0) write (*,*) '32 ',i,xhelp_a
+          xhelp = tw32(left,i)/(tof32(left,i,iadc)**0.5)
+          tof32(left,i,itdc) = tof32(left,i,itdc) + xhelp
-        if(xhelp_a<3786) xhelp = tw32(left,i)/sqrt(xhelp_a)
+          tdc_c(ch32a(i),hb32a(i))=tof32(left,i,itdc)
-        tof32(left,i,itdc) = xhelp_t  + xhelp
+                                               ENDIF
-        tdc_c(ch32a(i),hb32a(i))=tof32(left,i,itdc)
-        xhelp_a = tof32(right,i,iadc)
+          if ((tdc(ch32b(i),hb32b(i)).lt.4095).and.
-        xhelp_t = tof32(right,i,itdc)
+      &             (tof32(right,i,iadc).lt.3786)) THEN
-        if(xhelp_a<3786) xhelp = tw32(right,i)/sqrt(xhelp_a)
+          xhelp = tw32(right,i)/(tof32(right,i,iadc)**0.5)
-        tof32(right,i,itdc) = xhelp_t  + xhelp
+          tof32(right,i,itdc) = tof32(right,i,itdc) + xhelp
-        tdc_c(ch32b(i),hb32b(i))=tof32(right,i,itdc)
+          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
  C----------------------------------------------------------------------
-Line 883 
 c       DATA tof31_x/ -6.0,0.,6.0/
+Line 1083 
 c       DATA tof31_x/ -6.0,0.,6.0/
  c       DATA tof32_y/ -5.0,0.0,5.0/
        yhelp=0.
+ c      yhelp=100.
        if (tof12_i.GT.none_find) yhelp=tof12_y(tof12_i)
        if (ytofpos(1).lt.100)  yhelp=ytofpos(1)
-Line 893 
 c       DATA tof32_y/ -5.0,0.0,5.0/
+Line 1094 
 c       DATA tof32_y/ -5.0,0.0,5.0/
  c          if (adc(ch11a(i),hb11a(i)).lt.4095) then
              tof11(left,i,iadc) = tof11(left,i,iadc)*cos(theta13)
              xkorr = atten(left,11,i,yhelp)
- c            write(40+i,*) yhelp,xkorr
              xkorr=xkorr/hepratio
              adctof_c(ch11a(i),hb11a(i))=tof11(left,i,iadc)/xkorr
           endif
-Line 902 
 c            write(40+i,*) yhelp,xkorr
+Line 1102 
 c            write(40+i,*) yhelp,xkorr
  c          if (adc(ch11b(i),hb11b(i)).lt.4095) then
              tof11(right,i,iadc) = tof11(right,i,iadc)*cos(theta13)
              xkorr = atten(right,11,i,yhelp)
- c            write(40+i,*) yhelp,xkorr
              xkorr=xkorr/hepratio
              adctof_c(ch11b(i),hb11b(i))=tof11(right,i,iadc)/xkorr
           endif
        ENDIF
        xhelp=0.
+ c      xhelp=100.
        if (tof11_i.GT.none_find) xhelp=tof11_x(tof11_i)
        if (xtofpos(1).lt.100)  xhelp=xtofpos(1)
-Line 919 
 c            write(40+i,*) yhelp,xkorr
+Line 1119 
 c            write(40+i,*) yhelp,xkorr
  c          if (adc(ch12a(i),hb12a(i)).lt.4095) then
              tof12(left,i,iadc) = tof12(left,i,iadc)*cos(theta13)
              xkorr = atten(left,12,i,xhelp)
- c            write(50+i,*) xhelp,xkorr
              xkorr=xkorr/hepratio
              adctof_c(ch12a(i),hb12a(i))=tof12(left,i,iadc)/xkorr
           endif
-Line 928 
 c            write(50+i,*) xhelp,xkorr
+Line 1127 
 c            write(50+i,*) xhelp,xkorr
  c          if (adc(ch12b(i),hb12b(i)).lt.4095) then
              tof12(right,i,iadc) = tof12(right,i,iadc)*cos(theta13)
              xkorr = atten(right,12,i,xhelp)
- c            write(50+i,*) xhelp,xkorr
              xkorr=xkorr/hepratio
              adctof_c(ch12b(i),hb12b(i))=tof12(right,i,iadc)/xkorr
           endif
-Line 937 
 c            write(50+i,*) xhelp,xkorr
+Line 1135 
 c            write(50+i,*) xhelp,xkorr
  C-----------------------------S2 --------------------------------
        xhelp=0.
+ c      xhelp=100.
        if (tof22_i.GT.none_find) xhelp=tof22_x(tof22_i)
        if (xtofpos(2).lt.100)  xhelp=xtofpos(2)
-Line 947 
 C-----------------------------S2 -------
+Line 1146 
 C-----------------------------S2 -------
  c          if (adc(ch21a(i),hb21a(i)).lt.4095) then
              tof21(left,i,iadc) = tof21(left,i,iadc)*cos(theta13)
              xkorr = atten(left,21,i,xhelp)
- c            write(60+i,*) xhelp,xkorr
              xkorr=xkorr/hepratio
              adctof_c(ch21a(i),hb21a(i))=tof21(left,i,iadc)/xkorr
           endif
-Line 957 
 c          if (adc(ch21b(i),hb21b(i)).lt
+Line 1155 
 c          if (adc(ch21b(i),hb21b(i)).lt
              tof21(right,i,iadc) = tof21(right,i,iadc)*cos(theta13)
              xkorr=adcx21(right,i,1)*exp(xhelp/adcx21(right,i,2))
              xkorr = atten(right,21,i,xhelp)
- c            write(60+i,*) xhelp,xkorr
              xkorr=xkorr/hepratio
              adctof_c(ch21b(i),hb21b(i))=tof21(right,i,iadc)/xkorr
           endif
-Line 965 
 c            write(60+i,*) xhelp,xkorr
+Line 1162 
 c            write(60+i,*) xhelp,xkorr
        yhelp=0.
+ c      yhelp=100.
        if (tof21_i.GT.none_find) yhelp=tof21_y(tof21_i)
        if (ytofpos(2).lt.100)  yhelp=ytofpos(2)
-Line 975 
 c            write(60+i,*) xhelp,xkorr
+Line 1173 
 c            write(60+i,*) xhelp,xkorr
  c          if (adc(ch22a(i),hb22a(i)).lt.4095) then
              tof22(left,i,iadc) = tof22(left,i,iadc)*cos(theta13)
              xkorr = atten(left,22,i,yhelp)
- c            write(70+i,*) yhelp,xkorr
              xkorr=xkorr/hepratio
              adctof_c(ch22a(i),hb22a(i))=tof22(left,i,iadc)/xkorr
           endif
-Line 984 
 c            write(70+i,*) yhelp,xkorr
+Line 1181 
 c            write(70+i,*) yhelp,xkorr
  c          if (adc(ch22b(i),hb22b(i)).lt.4095) then
              tof22(right,i,iadc) = tof22(right,i,iadc)*cos(theta13)
              xkorr = atten(right,22,i,yhelp)
- c            write(70+i,*) yhelp,xkorr
              xkorr=xkorr/hepratio
              adctof_c(ch22b(i),hb22b(i))=tof22(right,i,iadc)/xkorr
           endif
-Line 993 
 c            write(70+i,*) yhelp,xkorr
+Line 1189 
 c            write(70+i,*) yhelp,xkorr
  C-----------------------------S3 --------------------------------
        yhelp=0.
+ c      yhelp=100.
        if (tof32_i.GT.none_find) yhelp=tof32_y(tof32_i)
        if (ytofpos(3).lt.100)  yhelp=ytofpos(3)
-Line 1003 
 C-----------------------------S3 -------
+Line 1200 
 C-----------------------------S3 -------
  c          if (adc(ch31a(i),hb31a(i)).lt.4095) then
              tof31(left,i,iadc) = tof31(left,i,iadc)*cos(theta13)
              xkorr = atten(left,31,i,yhelp)
- c            write(80+i,*) yhelp,xkorr
              xkorr=xkorr/hepratio
              adctof_c(ch31a(i),hb31a(i))=tof31(left,i,iadc)/xkorr
           endif
-Line 1012 
 c            write(80+i,*) yhelp,xkorr
+Line 1208 
 c            write(80+i,*) yhelp,xkorr
  c          if (adc(ch31b(i),hb31b(i)).lt.4095) then
              tof31(right,i,iadc) = tof31(right,i,iadc)*cos(theta13)
              xkorr = atten(right,31,i,yhelp)
- c            write(80+i,*) yhelp,xkorr
              xkorr=xkorr/hepratio
              adctof_c(ch31b(i),hb31b(i))=tof31(right,i,iadc)/xkorr
           endif
        ENDIF
        xhelp=0.
+ c      xhelp=100.
        if (tof31_i.GT.none_find) xhelp=tof31_x(tof31_i)
        if (xtofpos(3).lt.100)  xhelp=xtofpos(3)
-Line 1029 
 c            write(80+i,*) yhelp,xkorr
+Line 1225 
 c            write(80+i,*) yhelp,xkorr
  c          if (adc(ch32a(i),hb32a(i)).lt.4095) then
              tof32(left,i,iadc) = tof32(left,i,iadc)*cos(theta13)
              xkorr = atten(left,32,i,xhelp)
- c            write(90+i,*) xhelp,xkorr
              xkorr=xkorr/hepratio
              adctof_c(ch32a(i),hb32a(i))=tof32(left,i,iadc)/xkorr
           endif
-Line 1038 
 c            write(90+i,*) xhelp,xkorr
+Line 1233 
 c            write(90+i,*) xhelp,xkorr
  c          if (adc(ch32b(i),hb32b(i)).lt.4095) then
              tof32(right,i,iadc) = tof32(right,i,iadc)*cos(theta13)
              xkorr = atten(right,32,i,xhelp)
- c            write(90+i,*) xhelp,xkorr
              xkorr=xkorr/hepratio
              adctof_c(ch32b(i),hb32b(i))=tof32(right,i,iadc)/xkorr
           endif
        ENDIF
  C--------------------------------------------------------------------
  C----------------------calculate Beta  ------------------------------
  C--------------------------------------------------------------------
-Line 1065 
 C     S11 - S31
+Line 1258 
 C     S11 - S31
           ds = xhelp1-xhelp2
           ihelp=(tof11_i-1)*3+tof31_i
           c1 = k_S11S31(1,ihelp)
+          if (iz.gt.2) c1 = c1 + k1corrA1
           c2 = k_S11S31(2,ihelp)
           betatof_a(1) = c2/(cos(theta13)*(ds-c1))
-Line 1093 
 C     S11 - S32
+Line 1287 
 C     S11 - S32
           ds = xhelp1-xhelp2
           ihelp=(tof11_i-1)*3+tof32_i
           c1 = k_S11S32(1,ihelp)
+          if (iz.gt.2) c1 = c1 + k1corrA1
           c2 = k_S11S32(2,ihelp)
           betatof_a(2) = c2/(cos(theta13)*(ds-c1))
-Line 1121 
 C     S12 - S31
+Line 1316 
 C     S12 - S31
           ds = xhelp1-xhelp2
           ihelp=(tof12_i-1)*3+tof31_i
           c1 = k_S12S31(1,ihelp)
+          if (iz.gt.2) c1 = c1 + k1corrA1
           c2 = k_S12S31(2,ihelp)
           betatof_a(3) = c2/(cos(theta13)*(ds-c1))
-Line 1149 
 C     S12 - S32
+Line 1345 
 C     S12 - S32
           ds = xhelp1-xhelp2
           ihelp=(tof12_i-1)*3+tof32_i
           c1 = k_S12S32(1,ihelp)
+          if (iz.gt.2) c1 = c1 + k1corrA1
           c2 = k_S12S32(2,ihelp)
           betatof_a(4) = c2/(cos(theta13)*(ds-c1))
-Line 1177 
 C     S21 - S31
+Line 1374 
 C     S21 - S31
           ds = xhelp1-xhelp2
           ihelp=(tof21_i-1)*3+tof31_i
           c1 = k_S21S31(1,ihelp)
+          if (iz.gt.2) c1 = c1 + k1corrB1
           c2 = k_S21S31(2,ihelp)
           betatof_a(5) = c2/(cos(theta13)*(ds-c1))
-Line 1205 
 C     S21 - S32
+Line 1403 
 C     S21 - S32
           ds = xhelp1-xhelp2
           ihelp=(tof21_i-1)*3+tof32_i
           c1 = k_S21S32(1,ihelp)
+          if (iz.gt.2) c1 = c1 + k1corrB1
           c2 = k_S21S32(2,ihelp)
           betatof_a(6) = c2/(cos(theta13)*(ds-c1))
-Line 1233 
 C     S22 - S31
+Line 1432 
 C     S22 - S31
           ds = xhelp1-xhelp2
           ihelp=(tof22_i-1)*3+tof31_i
           c1 = k_S22S31(1,ihelp)
+          if (iz.gt.2) c1 = c1 + k1corrB1
           c2 = k_S22S31(2,ihelp)
           betatof_a(7) = c2/(cos(theta13)*(ds-c1))
-Line 1261 
 C     S22 - S32
+Line 1461 
 C     S22 - S32
           ds = xhelp1-xhelp2
           ihelp=(tof22_i-1)*3+tof32_i
           c1 = k_S22S32(1,ihelp)
+          if (iz.gt.2) c1 = c1 + k1corrB1
           c2 = k_S22S32(2,ihelp)
           betatof_a(8) = c2/(cos(theta13)*(ds-c1))
-Line 1289 
 C     S11 - S21
+Line 1490 
 C     S11 - S21
           ds = xhelp1-xhelp2
           ihelp=(tof11_i-1)*2+tof21_i
           c1 = k_S11S21(1,ihelp)
+          if (iz.gt.2) c1 = c1 + k1corrC1
           c2 = k_S11S21(2,ihelp)
           betatof_a(9) = c2/(cos(theta13)*(ds-c1))
-Line 1317 
 C     S11 - S22
+Line 1519 
 C     S11 - S22
           ds = xhelp1-xhelp2
           ihelp=(tof11_i-1)*2+tof22_i
           c1 = k_S11S22(1,ihelp)
+          if (iz.gt.2) c1 = c1 + k1corrC1
           c2 = k_S11S22(2,ihelp)
           betatof_a(10) = c2/(cos(theta13)*(ds-c1))
-Line 1345 
 C     S12 - S21
+Line 1548 
 C     S12 - S21
           ds = xhelp1-xhelp2
           ihelp=(tof12_i-1)*2+tof21_i
           c1 = k_S12S21(1,ihelp)
+          if (iz.gt.2) c1 = c1 + k1corrC1
           c2 = k_S12S21(2,ihelp)
           betatof_a(11) = c2/(cos(theta13)*(ds-c1))
-Line 1373 
 C     S12 - S22
+Line 1577 
 C     S12 - S22
           ds = xhelp1-xhelp2
           ihelp=(tof12_i-1)*2+tof22_i
           c1 = k_S12S22(1,ihelp)
+          if (iz.gt.2) c1 = c1 + k1corrC1
           c2 = k_S12S22(2,ihelp)
           betatof_a(12) = c2/(cos(theta13)*(ds-c1))
-Line 1393 
 C-------
+Line 1598 
 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
-       icount=0
+ C--------  New mean beta  calculation  -----------------------
-       sw=0.
-       sxw=0.
-       beta_mean=100.
-       do i=1,12
+         do i=1,12
-          if ((betatof_a(i).gt.-1.5).and.(betatof_a(i).lt.1.5)) then
+          btemp(i) = betatof_a(i)
-             icount= icount+1
+         enddo
-             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 + betatof_a(i)*w_i
-             sw =sw + w_i
-          endif
-       enddo
-       if (icount.gt.0) beta_mean=sxw/sw
-       betatof_a(13) = beta_mean
+         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
-Line 1423 
 c      write(*,*) ytofpos
+Line 1639 
 c      write(*,*) ytofpos
  c      write(*,*) xtr_tof
  c      write(*,*) ytr_tof
- continue
+ c 100  continue
+         continue
  C
        RETURN
-Line 1515 
 c       write(*,*) ix,pc_adc
+Line 1732 
 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: ',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
+         REAL tdcfl(4,12)
+         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****************************************************************************

 Legend:



Removed from v.1.6
 


changed lines


 
Added in v.1.12
 Legend:



Removed from v.1.6
 


changed lines


 
Added in v.1.12
-Removed from v.1.6
+Added in v.1.12

	ViewVC Help
Powered by ViewVC 1.1.23