MODULE MYanalyzer
   USE IOunits
   USE Threshold
   USE TPCdata
   
   USE TPCpaw
   USE TPCcalcAmplitude
   USE TPCtrack
   USE MyMinuit
   USE Space
   
   PRIVATE 
   
   PUBLIC  :: DenseData
   PUBLIC  :: analyzEvent,MY_book
   PRIVATE :: NF_analysis,HF_test,HF_over
   PRIVATE :: noOverflow,determinePRF,fitAmplRow,PRFCHI
   
CONTAINS
   
   SUBROUTINE analyzEvent(iEvent)
      INTEGER,INTENT(IN) :: iEvent
      INTEGER :: IERR,IVETO,nOver,iGmax
      REAL,DIMENSION(100) :: Xdata
      TYPE(Line_type) :: seedTrack,ZTrack,fitTrack,rowTrack
      REAL(KIND=spPrec) :: Z,THETA,X0Start,PhiStart,X0Fit,PhiFit,Width
      REAL :: dX2,dX3,dX4,dX5,dX6,dXO3,dXO4,dXO5,T0,T02,T03,T04,T05,T06, &
              AmplRow2,AmplRow3,AmplRow4,AmplRow5,AmplRow6
      INTEGER,DIMENSION(10) :: iGroupV=(/8,14,20,26,32,34,40,46,52,58/)
      INTEGER,DIMENSION(1) :: iRowSeed1=(/2/),iRowSeed2=(/6/)
      INTEGER,DIMENSION(5) :: iRowT=(/2,3,4,5,6/)
      INTEGER,DIMENSION(4) :: iRowTT
      INTEGER,DIMENSION(1) :: iRowF

! Where am I?
!      PRINT*,"Event: ",iEvent,EventNumber,"  Run:",RunNumber
     
! get the signal amplitudes
      CALL calcAmplitudeFoil() 
!      CALL HF_amplitude(iRowT)

! preselection
      CALL eventVeto(iGroupV,iRowSeed1,iRowSeed2,iRowT,IVETO)
      CALL HF_veto(IVETO)
      IF( IVETO /= 0 )RETURN

      nOver = noOverflow(iRowT)
      CALL HF_over(nOver)
      IF( nOver /= 0 )RETURN

! correct the gain
!      CALL correctGain()  

! plain event dispay and plot pulse shapes
!      CALL NF_event()
!      CALL plotPulseRow(iRowT)
!      CALL plotPulseEvent()
!      call hf_test(iEvent)

! find seed track and NF for event display
      CALL findSeedTrack(iRowSeed1,iRowSeed2,iRowT,seedTrack,IERR)
      CALL HF_seedTrack(IERR+1)
   !  CALL NF_event(seedTrack)

      IF( IERR /= 0 ) RETURN

! Z-Y track, linear regression or fit
!      CALL TimeLine(iRowT,time0,deltaT)
      CALL TrackFitZ(iRowT,ZTrack)
      CALL getLine(ZTrack,Z,THETA)
      Width = trackWidth(ZTrack,0.0)

      fitTrack = seedTrack
! X-Y track, with or without variable track width, NF for event display
!     fit 3 parameters
!      CALL TrackFitWidth(iRowT,fitTrack,Width)
!     fit 2 parameters, calculate width from time0
!      CALL TrackFitXY(iRowT,time0,fitTrack)
      CALL TrackFitXY(iRowT,ZTrack,fitTrack)
      CALL NF_event(fitTrack,Width)
      CALL getLine(seedTrack,X0Start,PhiStart)
      CALL getLine(fitTrack,X0Fit,PhiFit)
      
      DX2 = 1000.0
      DX3 = 1000.0
      DX4 = 1000.0
      DX5 = 1000.0
      DX6 = 1000.0
      DXO3 = 1000.0
      DXO4 = 1000.0
      DXO5 = 1000.0
      
      AmplRow3 = 0.0
      AmplRow4 = 0.0
      AmplRow5 = 0.0
      
      CALL GmaxSignalRow(2,iGmax,AmplRow2,T02)
      CALL GmaxSignalRow(3,iGmax,AmplRow3,T03)
      CALL GmaxSignalRow(4,iGmax,AmplRow4,T04)
      CALL GmaxSignalRow(5,iGmax,AmplRow5,T05)
      CALL GmaxSignalRow(6,iGmax,AmplRow6,T06)
      
      
      IF( .NOT.(fitTrack == seedTrack) )THEN
! resolution fit
         iRowF(1) = 2
         CALL TrackFitX(iRowF,Width,fitTrack,dX2)
         iRowF(1) = 3
         CALL TrackFitX(iRowF,Width,fitTrack,dX3)
         iRowF(1) = 4
         CALL TrackFitX(iRowF,Width,fitTrack,dX4)
         iRowF(1) = 5
         CALL TrackFitX(iRowF,Width,fitTrack,dX5)
         iRowF(1) = 6
         CALL TrackFitX(iRowF,Width,fitTrack,dX6)
         
         rowTrack = fitTrack
         iRowTT = (/2,4,5,6/)
         iRowF(1) = 3
         CALL TrackFitXY(iRowTT,ZTrack,rowTrack)
         CALL TrackFitX(iRowF,Width,rowTrack,dXO3)
         rowTrack = fitTrack
         iRowTT = (/2,3,5,6/)
         iRowF(1) = 4
         CALL TrackFitXY(iRowTT,ZTrack,rowTrack)
         CALL TrackFitX(iRowF,Width,rowTrack,dXO4)
         rowTrack = fitTrack
         iRowTT = (/2,3,4,6/)
         iRowF(1) = 5
         CALL TrackFitXY(iRowTT,ZTrack,rowTrack)
         CALL TrackFitX(iRowF,Width,rowTrack,dXO5)

         CALL fitAmplRow(fitTrack,Z,3,AmplRow3)
         CALL fitAmplRow(fitTrack,Z,4,AmplRow4)
         CALL fitAmplRow(fitTrack,Z,5,AmplRow5)
      ELSE
! fit failed      
         X0Fit = 10000.0
         PhiFit = 10000.0
      ENDIF
      
! fill ntuple for analysis   

      Xdata = -1.0E6
      Xdata(1) = RunNumber
      Xdata(2) = EventNumber
      Xdata(3) = X0Start
      Xdata(4) = PhiStart
      Xdata(5) = X0Fit
      Xdata(6) = PhiFit
      Xdata(7) = Width
      Xdata(8) = Z
      Xdata(9) = Theta
      Xdata(10) = dXO3
      Xdata(11) = dXO4
      Xdata(12) = dXO5
      Xdata(13) = dX2
      Xdata(14) = dX3
      Xdata(15) = dX4
      Xdata(16) = dX5
      Xdata(17) = dX6
      Xdata(18) = AmplRow3
      Xdata(19) = AmplRow4
      Xdata(20) = AmplRow5
      CALL NF_analysis(Xdata)
      
!      CALL determinePRF(fitTrack,Z)
   END SUBROUTINE analyzEvent
   
   FUNCTION noOverflow(iRowT) RESULT(nOver)
      INTEGER,DIMENSION(:),INTENT(IN) :: iRowT
      INTEGER :: nOver,iR,iRow,iG,iGroup
   
      nOver = 0
! reject events with overflow
      DO iR=1,size(iRowT)
         iRow = iRowT(iR)
         DO iG=1,nGroupsinRow(iRow)
            iGroup = iGroupofRow(iG,iRow)
            IF( minval(IADC(1:nTbGroup(iGroup),iGroup)) <= 0 ) nOver = nOver+1
         ENDDO
      ENDDO
   END FUNCTION noOverflow

   SUBROUTINE determinePRF(fitTrack,Z)
      TYPE(Line_type),INTENT(IN) :: fitTrack
      REAL(KIND=spPrec),INTENT(IN) :: Z
      INTEGER :: iR,iRow,iP,iPad,iGmax
      REAL :: Amplitude,amplrow,Xpad,Ypad,xTrack,deltaX,deltaT
      REAL :: T0,Trow,AmplMax
      REAL,DIMENSION(100) :: Xdata
      INTEGER,DIMENSION(3) :: iRowTest=(/3,4,5/)
      TYPE(Coordinate_type) :: Point     
      REAL(KIND=spPrec) :: X0Fit,PhiFit
     
      CALL getLine(fitTrack,X0Fit,PhiFit)
     
      DO iR = 1,size(iRowTest)
         iRow = iRowTest(iR)
         CALL GmaxSignalRow(iRow,iGmax,AmplMax,Trow)
         IF( AmplMax>200.0 ) CYCLE

         CALL GetSignalRow(iRow,AmplRow,T0)
         IF( AmplRow<30.0 ) CYCLE
    
         CALL fitAmplRow(fitTrack,Z,iRow,AmplRow) 
    
         DO iP = 1,nPadsinRow(iRow)
            iPad = iPadofRow(iP,iRow)
            CALL GetSignal(iPad,Amplitude,T0)

            Point = getCentre(PadLocation(iPad))
            Xpad = Pitch(Point)
            Ypad = Length(Point)
            Xtrack = getLineX(fitTrack,Ypad)
            deltaX = Distance(fitTrack,Point)
            IF( Xpad<Xtrack ) deltaX = -deltaX

!                DD = getLineX(ZTrack,getCentreY(iPad))
!                tTrack = DD / tBin2ns / vDrift * 0.1 + tBin0
!                deltaT = T0 - tTrack
            deltaT = T0 - Trow
            
            IF( Amplitude > 0.001 .OR. ABS(deltaX)<6.0 )THEN 
               XDATA(1) = X0Fit
               XDATA(2) = PhiFit
               XDATA(3) = Z
               XDATA(4) = deltaX
               XDATA(5) = deltaT
               XDATA(6) = Amplitude/AmplRow
               XDATA(7) = AmplRow
               XDATA(8) = iPad
               CALL HFN(952,Xdata)
            ENDIF      
         ENDDO
      ENDDO

   END SUBROUTINE determinePRF
   SUBROUTINE fitAmplRow(fitTrack,Z,iRow,AmplRow)
      TYPE(Line_type),INTENT(IN) :: fitTrack
      REAL(KIND=spPrec),INTENT(IN) :: Z
      INTEGER,INTENT(IN) :: iRow
      REAL,INTENT(OUT) :: AmplRow
      INTEGER :: IPRT,IERR
      REAL(KIND=spPrec) :: X0,Phi
      CHARACTER(LEN=50) :: CHTITL
      INTEGER,PARAMETER :: NFITP=5
      CHARACTER(LEN=10),DIMENSION(NFITP) :: PCH
      REAL(KIND=minuitPrec) :: CHISQ
      REAL(KIND=minuitPrec),DIMENSION(NFITP) :: PVAL,PERR,BND1,BND2
      REAL(KIND=minuitPrec),DIMENSION(NFITP,NFITP) :: PMAT
     
      CALL getLine(fitTrack,X0,Phi)
     
      PCH(1) = "AmplRow   "
      PCH(2) = "X0        "
      PCH(3) = "Phi       "
      PCH(4) = "Sigma     "
      PCH(5) = "Row       "
      PVAL(1) = 30.0      
      PVAL(2) = X0
      PVAL(3) = Phi
      PVAL(4) = trackWidth(Line(Z,0.0_spPrec),0.0)     
      PVAL(5) = iRow      
      PERR(1) = 3.0
      PERR(2) = -0.1
      PERR(3) = -1.0
      PERR(4) = -1.0
      PERR(5) = -1.0
! normaly no bounds
      BND1 = 0.0      
      BND2 = 0.0  
      WRITE(UNIT=CHTITL,FMT=*) "+++++ PRF fit, row:",iRow
      IPRT = 0
      CALL MINUIT(PRFCHI,CHTITL,NFITP,PCH,PVAL,PERR,PMAT, &
                     BND1,BND2,CHISQ,IPRT,IERR)
      AmplRow = PVAL(1) 
    
   END SUBROUTINE fitAmplRow
   
   SUBROUTINE PRFCHI(NPAR,GRAD,FVAL,XVAL,IFLAG)
      INTEGER,INTENT(IN) :: IFLAG , NPAR
      REAL(KIND=minuitPrec),DIMENSION(5),INTENT(IN) :: XVAL
      REAL(KIND=minuitPrec),DIMENSION(5),INTENT(OUT) :: GRAD
      REAL(KIND=minuitPrec),INTENT(OUT) :: FVAL
      REAL(KIND=minuitPrec) :: CHISQ,sum
      TYPE(Line_type) :: Track   
      INTEGER :: iRow,iP,iPad,myPar
      REAL(KIND=spPrec) :: X0,Phi,width
      REAL :: AmplRow,Amplitude,T0

! dummy to avoid compiler warnings
      IF( IFLAG==9 ) myPar=NPAR
      
      GRAD   = 0.0_minuitPrec
      CHISQ = 0.0_minuitPrec  
      
      AmplRow = XVAL(1)
      X0      = XVAL(2)
      Phi     = XVAL(3)
      width   = XVAL(4)
      iRow    = XVAL(5)
      Track   = Line(X0,Phi)
         
      DO iP = 1,nPadsinRow(iRow)
         iPad = iPadofRow(iP,iRow)
         CALL GetSignal(iPad,Amplitude,T0)
         
         sum = (Amplitude-PRF(Track,width,iPad)*AmplRow)/AERR(Amplitude)
         CHISQ = CHISQ + sum**2
      ENDDO
      
      FVAL = CHISQ
   END SUBROUTINE PRFCHI

   SUBROUTINE MY_book()
! book histograms & Ntuples for analyzEvent
      CHARACTER(LEN=8),DIMENSION(1000) :: CHtag
      CHARACTER(LEN=20) :: CHtitle
      INTEGER :: ntag,i
      
! Histograms, numbers >1039  
      CALL HBOOK1(1040,"test",64,0.5,64.5,0.0)     
      call hbook1(1041,"t0",400,0.0,4000.0,0.0)
      call hbook1(1050,"nOver",40,-0.5,39.5,0.0)
      
!      do i=1,ngroups
!         WRITE(UNIT=CHtitle,FMT=*) "mean amplitude",i
!         CALL HBPROF(1100+i,CHtitle,50,0.0,10000.0,-0.5,4000.0,' ')
!      enddo
      
      CALL HBOOK2(1046,"time difference",40,0.0,0.5,40,-50.0,350.0,0.0)
      CALL HBOOK2(1047,"shape",40,0.0,0.5,40,0.0,1.0,0.0)
      CALL HBOOK2(1048,"shape",40,-50.0,350.0,40,0.0,1.0,0.0)

      CHtag(1)  = "RUN"
      CHtag(2)  = "EVENT"
      CHtag(3)  = "X0START"
      CHtag(4)  = "PHISTART"
      CHtag(5)  = "X0FIT"
      CHtag(6)  = "PHIFIT"
      CHtag(7)  = "SIGMAFIT"
      CHtag(8)  = "Z"
      CHtag(9)  = "THETA"
      CHtag(10) = "DXWOUT3"
      CHtag(11) = "DXWOUT4"
      CHtag(12) = "DXWOUT5"
      CHtag(13) = "DXWITH2"
      CHtag(14) = "DXWITH3"
      CHtag(15) = "DXWITH4"
      CHtag(16) = "DXWITH5"
      CHtag(17) = "DXWITH6"
      CHtag(18) = "AMPLR3"
      CHtag(19) = "AMPLR4"
      CHtag(20) = "AMPLR5"
      ntag = 20
      CALL NT_book(951,"Analysis",ntag,CHtag)
  
      CHtag(1) = "X0FIT"
      CHtag(2) = "PHIFIT"
      CHtag(3) = "Z"
      CHtag(4) = "DELTAX"
      CHtag(5) = "DELTAT"
      CHtag(6) = "AMPLNORM"
      CHtag(7) = "AMPLROW"
      CHtag(8) = "PAD"
      ntag = 8
      CALL NT_book(952,"PRF",ntag,CHtag)
  

   END SUBROUTINE MY_book
! NTUPLE fill  ********************************************** 
   SUBROUTINE NF_analysis(Xdata)
      REAL,DIMENSION(:),INTENT(IN) :: Xdata
      
      CALL HFN(951,Xdata)
   END SUBROUTINE NF_analysis
      
! HISTOGRAM fill  ********************************************** 
   SUBROUTINE HF_test(iEvent)
      INTEGER,INTENT(IN) :: iEvent
      integer:: iRow,iG,iGroup,nHitRow,iTmax,iT0
      real :: amplitude,t0,amplsum,amplratio,ampllate,amplmax
      
      DO iGroup=1,nGroups
         CALL getSignalGroup(iGroup,Amplitude,T0)      
         if( amplitude>1.0 )then
         CALL HFILL(1040,real(iGroup),0.0,amplitude)
         call hfill(1041,t0,0.0,1.0)
         endif
!         call hfill(1100+iGroup,REAL(iEvent),amplitude,1.0)
      enddo

      DO iRow=1,nRows
         amplsum = 0.0
         amplmax = 0.0
         nHitRow = 0
         DO iG=1,nGroupsinRow(iRow)
            iGroup = iGroupofRow(iG,iRow)
            CALL getSignalGroup(iGroup,Amplitude,T0)      
            amplsum = amplsum + Amplitude
            IF( Amplitude > amplmax )THEN
               iTmax = T0
               amplmax = Amplitude
            ENDIF
            IF( Amplitude > 0.0 ) nHitRow = nHitRow + 1
         ENDDO
         IF( nHitRow > 1 )THEN
         DO iG=1,nGroupsinRow(iRow)
            iGroup = iGroupofRow(iG,iRow)
            CALL getSignalGroup(iGroup,Amplitude,T0)   
            IF( Amplitude > 0.0 )THEN
            amplratio = Amplitude / amplsum
            
            iT0 = T0
            IF( nTbin/=0 )THEN
            ampllate = sum(ADC(iT0+81:iT0+130,iGroup)) / 50.0
            ampllate = -ampllate / Amplitude
            ELSE
               ampllate = 0.0
            ENDIF
            
            CALL HF2(1046,amplratio,REAL(iT0-iTmax),1.0)
            CALL HF2(1047,amplratio,ampllate,1.0)
            CALL HF2(1048,REAL(iT0-iTmax),ampllate,1.0)
            ENDIF
         ENDDO
         ENDIF
      ENDDO 

   END SUBROUTINE HF_test  
   SUBROUTINE HF_over(nOver)
      INTEGER,INTENT(IN) :: nOver
      
      CALL HF1(1050,REAL(nOver),1.0)
   END SUBROUTINE HF_over

! >>>>  WRITE DENSE DATA  >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>

   SUBROUTINE DenseData()
   USE DataFormats, only : WriteDenseData
!      INTEGER :: IVETO

! get the signal amplitudes
      CALL calcAmplitudeFoil() 

! preselection
!      CALL eventVeto(IVETO)
!      IF( IVETO /= 0 )RETURN

      CALL WriteDenseData()
      
   END SUBROUTINE DenseData

END MODULE MYanalyzer