C------------------------------------------------
      PROGRAM TOFTRACK
C------------------------------------------------
C     W. Menn
C
C     Based on the "template" progran from  Elena Vannuccini
C
C     Version 1.00  August 2005
C     Version 1.01  07-oct-2005: Some bugs found
C     - new rz file must be named "***_toftrack.rz'
C     - booking of the EVENT and CPU blocks changed according
C       to Elena's "book_level2" routine
C      
C------------------------------------------------

      include 'trk_level2.f'
      include 'common_tofroutine.f'
      include 'common_tof.f'

c----- HBOOK
      INTEGER HMEM
      parameter (NWPAWC=8500000)
      common/PAWC/HMEM(NWPAWC)

      parameter (ntp_level2=22)
      parameter (ntp_new=2)

c------------------------------------------------------------------------
c     local variables
c------------------------------------------------------------------------
      character*74 data_file    !data file name
      character*74 tof_file     !tofdata file name
      character*74 data_dir     !data directory
      character*74 data_file_level2
      character*74 data_file_tof
      character*74 data_file_tofcalib
      character*74 data_file_new

      parameter (lun_data_level2=72) !data file id number

c------------------------------------------------------------------
      COMMON/QUEST/IQUEST(100)

C------------------------------------------------------------------
c     =======================================
c     variables for tracking routine
c     =======================================
      parameter(NPOINT_MAX=100)
      DOUBLE PRECISION ZIN(NPOINT_MAX)
      DOUBLE PRECISION XOUT1(NPOINT_MAX),YOUT1(NPOINT_MAX)
      REAL  XOUT(NPOINT_MAX),YOUT(NPOINT_MAX)
      DOUBLE PRECISION AL_P(5)
      REAL XFIT(6),YFIT(6),ZFIT(6),ALPHA(5),XGOOD1(6),YGOOD1(6)
c     =======================================

c     define TOF Z-coordinates
      parameter (NPTOF=3)
      REAL ZTOF(NPTOF)
c      DATA ZTOF/55.,25.,-23/    !NB!!!  from TOP to BOTTOM
      DATA ZTOF/53.4,23.9,-23.4/  !Stockholm 29.9.2005

c------------------------------------------------------------------------
C
      CHARACTER*100 name
C

      INTEGER lrec, openntuple
      INTEGER istat, ierr, icycle
      PARAMETER (lrec=4096)


      integer patterntrig1(6)
      real chi

      character*35 block1,block2

C-------------------------------------------------------------------------
c
c     HBOOK initialization
c
c------------------------------------------------------------------------

      call HLIMIT(NWPAWC)

C     HERE A PATH where you can write!
C
c      print *
c      print *,'File saved in: ',name
C
C     CHANGE here (hbname and hbnt if you want another ntuple id number
C

111    format(a)
       print*,'data dir:'
      read(*,111)data_dir
c       read(*,*)data_dir
      print*,data_dir
      print*,'tof file?'
      print*,'(without estention: output_YYMMDD_XXX )'
      read(*,*)tof_file

      print*,tof_file
      print*,' '

      print*,'first event:'
      read(*,*) minevent
      print*,minevent
c     print*,' '

      print*,'number of events to be analysed'
      read(*,*) ntotev
      print*,ntotev


c------------------------------------------------------------------------
c---------------    open ToF rz file 
c------------------------------------------------------------------------

      print*,
     $     data_dir(1:LNBLNK(data_dir))//tof_file(1:LNBLNK(tof_file))


505   format(a,'DW_',a,'_tof.rz')
      write(data_file_tof,505)
     $     data_dir(1:LNBLNK(data_dir))
     $     ,tof_file(1:LNBLNK(tof_file))
      print*,'__________ opening TOF rz file __________'
      print*,data_file_tof
 

      call HROPEN(59,
     $     'EVENT',data_file_tof,'QP',4096,istat) !opens rz


      if(istat.ne.0) stop
      print*,' '
      print*,'reading TOF n-tuple...'
      call HRIN(ntp_tof,9999,0)
      CALL HBNAME(ntp_tof,' ',0,'$CLEAR')
      CALL HBNAME(ntp_tof,'EVENT',good,'$SET')
      CALL HBNAME(ntp_tof,'TRIGGER',trig_evcount,'$SET')
      CALL HBNAME(ntp_tof,'TOF',tdcid,'$SET')
      call HPRNTU(ntp_tof)
      call HNOENT(ntp_tof,iemax0)
      write(*,*) 'Number of events ToF ',iemax0
 

c------------------------------------------------------------------------
c---------------    open ToF calib rz file 
c------------------------------------------------------------------------

      print*,
     $     data_dir(1:LNBLNK(data_dir))//tof_file(1:LNBLNK(tof_file))

506   format(a,'DW_',a,'_tofcalib.rz')
      write(data_file_tofcalib,506)
     $     data_dir(1:LNBLNK(data_dir))
     $     ,tof_file(1:LNBLNK(tof_file))
      print*,'__________ opening TOF calib rz file __________'
      print*,data_file_tofcalib
 
       CALL HROPEN(59,'TOF K1A',data_file_tofcalib,'QP',4096,istat)

        if (istat.ne.0) then  ! check if HROPEN was OK
              write(*,*) 'Can''t open correct ToF calibration File !!!'
              write(*,*) 'Will use standard file tofcalib.rz !!!'
507   format(a,'tofcalib.rz')
      write(data_file_tofcalib,507)
     $     data_dir(1:LNBLNK(data_dir))
      print*,'__________ opening TOF calib rz file __________'
      print*,data_file_tofcalib
 
       CALL HROPEN(59,'TOF K1A',data_file_tofcalib,'QP',4096,istat)
        endif

        print*,' reading TOF CALIB n-tuple...'

        call HRIN(ntp_tofcalib,9999,0)

        call HBNAME(ntp_tofcalib,' ',0,'$CLEAR')

        call HBNAME(ntp_tofcalib,'TOFK1A',k1_s11s31,'$SET')
        call HBNAME(ntp_tofcalib,'TOFK1B',k1_s12s32,'$SET')
        call HBNAME(ntp_tofcalib,'TOFK1C',k1_s21s31,'$SET')
        call HBNAME(ntp_tofcalib,'TOFK1D',k1_s22s32,'$SET')

        call HBNAME(ntp_tofcalib,'TOFLIN11',y_coor_lin11,'$SET')
        call HBNAME(ntp_tofcalib,'TOFLIN12',x_coor_lin12,'$SET')
        call HBNAME(ntp_tofcalib,'TOFLIN21',x_coor_lin21,'$SET')
        call HBNAME(ntp_tofcalib,'TOFLIN22',y_coor_lin22,'$SET')
        call HBNAME(ntp_tofcalib,'TOFLIN31',y_coor_lin31,'$SET')
        call HBNAME(ntp_tofcalib,'TOFLIN32',x_coor_lin32,'$SET')

        call HBNAME(ntp_tofcalib,'TOFTW11',tw11,'$SET')
        call HBNAME(ntp_tofcalib,'TOFTW12',tw12,'$SET')
        call HBNAME(ntp_tofcalib,'TOFTW21',tw21,'$SET')
        call HBNAME(ntp_tofcalib,'TOFTW22',tw22,'$SET')
        call HBNAME(ntp_tofcalib,'TOFTW31',tw31,'$SET')
        call HBNAME(ntp_tofcalib,'TOFTW32',tw32,'$SET')

        call HBNAME(ntp_tofcalib,'TOFADC11',adcx11,'$SET')
        call HBNAME(ntp_tofcalib,'TOFADC12',adcx12,'$SET')
        call HBNAME(ntp_tofcalib,'TOFADC21',adcx21,'$SET')
        call HBNAME(ntp_tofcalib,'TOFADC22',adcx22,'$SET')
        call HBNAME(ntp_tofcalib,'TOFADC31',adcx31,'$SET')
        call HBNAME(ntp_tofcalib,'TOFADC32',adcx32,'$SET')

        call HPRNTU(ntp_tofcalib)

        call HNOENT(ntp_tofcalib,iemax_cal)
c        write(*,*) 'Number of Events CALIB ',iemax_cal

        do iev=1,iemax_cal
        call HGNT(ntp_tofcalib,iev,ierr) !reads an event
        enddo

        call hrout(ntp_tofcalib,icycle,' ')
        call hrend('TOF K1A')
        close(0)



c------------------------------------------------------------------------
c---------------    open Track Level2 rz file 
c------------------------------------------------------------------------

      data_file=tof_file
      print*,'data dir:'
c      read(*,111)data_dir
      print*,data_dir
c      print*,'data file?'
c      print*,'(without estention: output_YYMMDD_XXX )'
c      read(*,*)data_file

      print*,'data_file:'
      print*,data_file
      print*,' '

      print*,
     $     data_dir(1:LNBLNK(data_dir))//data_file(1:LNBLNK(data_file))

c------------------------------------------------------------------------
 504  format(a,'DW_',a,'_level2.rz')
      write(data_file_level2,504)
     $     data_dir(1:LNBLNK(data_dir))
     $     ,data_file(1:LNBLNK(data_file))
      print*,'__________ opening LEVEL2 rz file __________'
      print*,data_file_level2
      print*,'__________ opening LEVEL2 rz file __________'
      print*,data_file_level2
      IQUEST(10)=65000
      call HROPEN(lun_data_level2,
     $     'LEVEL2',data_file_level2,'QP',4096,istat) !opens rz
      if(istat.ne.0) goto 19
      print*,'reading LEVEL2 n-tuple...'
      call HRIN(ntp_level2,9999,0)

*     -----------------------------------------------
      CALL HBNAME(ntp_level2,' ',0,'$CLEAR')
      CALL HBNAME(ntp_level2,'EVENT',GOOD2,'$SET')
      CALL HBNAME(ntp_level2,'CPU',pkt_type,'$SET')
      CALL HBNAME(ntp_level2,'TRACKS',ntrk,'$SET')
      CALL HBNAME(ntp_level2,'SINGLETS',nclsx,'$SET')
*     -----------------------------------------------

      call HPRNTU(ntp_level2)
      call HNOENT(ntp_level2,iemax2)

      write(*,*) 'Number of events LEVEL2 ',iemax2
      print*,'ok'
      print*,' '



c------------------------------------------------------------------------
c     read magnetic field map
c------------------------------------------------------------------------
c
c     =======================


c     read magnetic field map
c     =======================
c
c------------------------------------------------------------------------
      print*,'- read magnetic field map'
      print*,' '
      call read_B
      print*,' '

c      print*,' '
c      print*,'- read magnetic field map'
c      print*,' '
c      call read_B_2maps
c


c------------------------------------------------------------------------
c     open new RZ file    
c------------------------------------------------------------------------

503    format(a,'DW_',a,'_toftrack.rz')
       write(data_file_new,503)
     $     data_dir(1:LNBLNK(data_dir))
     $     ,data_file(1:LNBLNK(data_file))
       print*,'__________ opening NEW rz file __________'
       print*,data_file_new


      call HROPEN(10,'TEST',data_file_new,'NP',4096,istat) !opens rz

      call HBNT(ntp_new,'TOF',' ')


c      call HBNAME(ntp_new,'EVENT',good,'GOOD:L,NEV_TRK:I')

      call HBNAME(ntp_new,'EVENT',good2,'GOOD2:L,NEV2:I')


c      call HBNAME(ntp_new,'CPU',pkt_type
c     $     ,'PKT_TYPE:I
c     $     ,PKT_NUM:I
c     $     ,OBT:I
c     $     ,WHICH_CALIB:I')

      call HBNAME(ntp_new,'CPU',pkt_type
     $     ,'PKT_TYPE:I::[0,50]
     $     ,PKT_NUM:I
     $     ,OBT:I
     $     ,WHICH_CALIB:I::[0,50]')



      call HBNAME(ntp_new,'TOF',tdcid
     $     ,'TDCID(12):I
     $     ,EVCOUNT(12):I
     $     ,TDCMASK(12):I
     $     ,ADC(4,12):I
     $     ,TDC(4,12):I
     $     ,TEMP1(12):I
     $     ,TEMP2(12):I')

      call HBNAME(ntp_new,'TOF',beta_a,'BETA(5)')
      call HBNAME(ntp_new,'TOF',xtofpos,'XTOF(3)')
      call HBNAME(ntp_new,'TOF',ytofpos,'YTOF(3)')
      call HBNAME(ntp_new,'TOF',adc_c,'ADC_C(4,12)')
      call HBNAME(ntp_new,'TOF',tof_i_flag,'IFLAG(6)')
      call HBNAME(ntp_new,'TOF',tof_j_flag,'JFLAG(6)')

      call HBNAME(ntp_new,'TOF',xout,'XOUT(3)')
      call HBNAME(ntp_new,'TOF',yout,'YOUT(3)')


      call HBNAME(ntp_new,'TRIGGER',trig_evcount
     $     ,'TRIG_EVCOUNT:I
     $     ,PMTPL(3):I
     $     ,TRIGRATE(6):I
     $     ,DLTIME(2):I
     $     ,S4CALCOUNT(2):I
     $     ,PMTCOUNT1(24):I
     $     ,PMTCOUNT2(24):I
     $     ,PATTERNBUSY(3):I
     $     ,PATTERNTRIG(6):I
     $     ,TRIGCONF:I')


 417  format('NTRK:I::[0,',I4,']')
 418  format(',IMAGE(NTRK):I::[0,',I4,']')
      write(block1,417)NTRKMAX
      write(block2,418)NTRKMAX
      call HBNAME(ntp_new,'TRACKS',NTRK,
     $     block1//
     $     block2//'
     $     ,XM(6,NTRK):R
     $     ,YM(6,NTRK):R
     $     ,ZM(6,NTRK):R
     $     ,RESX(6,NTRK):R
     $     ,RESY(6,NTRK):R
     $     ,AL(5,NTRK):R
     $     ,COVAL(5,5,NTRK):R
     $     ,CHI2(NTRK):R
     $     ,XGOOD(6,NTRK):I::[0,1]
     $     ,YGOOD(6,NTRK):I::[0,1]
     $     ,XV(6,NTRK):R
     $     ,YV(6,NTRK):R
     $     ,ZV(6,NTRK):R
     $     ,AXV(6,NTRK):R
     $     ,AYV(6,NTRK):R
     $     ,DEDXP(6,NTRK):R
     $     ')
      call HBNAME(ntp_new,'SINGLETS',nclsx,
     $     'NCLSX(6):I,NCLSY(6):I')


      call HPRNTU(ntp_new)

c------------------------------------------------------------------------

      call HCDIR('//LEVEL2',' ')

c------------------------------------------------------------------------
c     start loop on events
c------------------------------------------------------------------------

       maxevent=minevent+ntotev
       igoodevent=0

c      do iev = 1,iemax0
       do iev = minevent,MIN(iemax0,maxevent) !loop on events      
       Call Hcdir('//EVENT',' ')
       Call Hgnt(ntp_tof,iev,Ierr)

               do i=1,6
                  patterntrig1(i) = patterntrig(i)
               enddo

*        ----------------------------------------------
       Call Hcdir('//LEVEL2',' ')
       call HGNT(ntp_level2,iev,ierr) !reads an event

       if(ierr.ne.0) goto 21
*        ----------------------------------------------

               do i=1,nptof
                  xout(i)=1000.
                  yout(i)=1000.
               enddo
               do i=1,5
                  ALPHA(i) = 1000.
               enddo
               do i=1,6
                  XFIT(i) = 1000.
                  YFIT(i) = 1000.
                  ZFIT(i) = 1000.
                  XGOOD1(i) = 1000.
                  YGOOD1(i) = 1000.
               enddo

               chi = 1000.

*=======> INSTERT HERE YOUR CODE

*        ---- example ----
*         write(*,*) 'GOOD2 ',good2
*         if(GOOD2)then
         if(NTRK.ne.0)then
            print*,'Event ',iev,' # tracks ',ntrk
            if(ntrk.eq.1)then

               chi = CHI2(1)  
               do i=1,5
                  ALPHA(i) = AL(i,1)
               enddo
               do i=1,6
                  XFIT(i) = XV(i,1)
                  YFIT(i) = YV(i,1)
                  ZFIT(i) = ZV(i,1)
                  XGOOD1(i)=XGOOD(i,1)
                  YGOOD1(i)=YGOOD(i,1)
               enddo

               igoodevent = igoodevent+1
*              assigned input  parameters for track routine
*              1) Z-coordinates where the trajectory is evaluated
               do itof=1,NPTOF
                  ZIN(itof) = ZTOF(itof)
               enddo
*              2) track status vector
               do i=1,5
                  AL_P(i) = AL(i,1)
               enddo
*              -------- *** tracking routine *** --------
               call track(NPTOF,ZIN,XOUT1,YOUT1,AL_P,IFAIL)
*              ------------------------------------------
               do itof=1,NPTOF
                  XOUT(itof)=XOUT1(itof)
                  YOUT(itof)=YOUT1(itof)
               enddo

               do itof=1,NPTOF
c                  print*,'   ',itof,ZIN(itof),XOUT(itof),YOUT(itof)
               enddo
            endif
         else
            print*,'Event ',iev,' *** BAD ***'
         endif


         call tofroutine(xout,yout,alpha)


c---------------------  fill new rz file   -----------------------------  

      call HFNT(2)

      enddo                     !end loop on events
      GOTO 9000                 !got to end

c------------------------------------------------------------------------
c
c     data file opening error
c
c------------------------------------------------------------------------
 19   continue

      print*,' '
      print*,'ERROR OPENING DATA FILE: ',data_file
      print*,' '
      print*,' '

      goto 9000                 !the end

c------------------------------------------------------------------------
c
c     level2 ntuple event reading error
c
c------------------------------------------------------------------------

 21   continue

      print*,' '
      print*,'ERROR WHILE READING LEVEL2 NTUPLE, AT EVENT
     $     : ',iev
      print*,' '
      print*,' '

      goto 9000                 !the end

c------------------------------------------------------------------------
c
c     closes files and exits
c
c------------------------------------------------------------------------

 9000 continue

      write(*,*) 'Good Events ',igoodevent
        

c--------------------------------------------------------------

      call HCDIR('//EVENT',' ')
      CALL HROUT(0,icycle,' ')
      call HREND('EVENT')
      close(0)
 
      call HCDIR('//LEVEL2',' ')
      call HREND('level2')
      close(lun_data_level2)


      call HPRNTU(2)
      call HCDIR('//TEST',' ')
      CALL HROUT(0,icycle,' ')
      call HREND('TEST')
      close(0)

      STOP
      END

 
      include 'tofroutine.f'

