MODULE TPCcalcAmplitude
   USE IOunits
   USE Threshold
   USE TPCdata
   USE TPCpaw
   USE MyMinuit  
   USE Space 
   PRIVATE 
   
! Calibration: determine pedestals, falltimes, ...
! Calibrate()
   PUBLIC :: Calibrate

! use data from TPCgain.txt to correct gain   
! correctGain()
   PUBLIC :: correctGain
   
! calcPedestal()
! subPedestal()
   PUBLIC :: calcPedestal,subPedestal

! ways to calculate the signal amplitude
! calcAmplitudeRow()
! calcAmplitudePad()
! calcAmplitudeFoil()
   PUBLIC :: calcAmplitudeRow,calcAmplitudePad,calcAmplitudeFoil
   PUBLIC :: calcAmplitudeFoilNEW

! to improve the amplitude for exponential pulse shapes
! fitAllPulses()
   PUBLIC :: fitAllPulses

! plot routines can use plotZoomPulse or plotFullPulse
! plot all pulses in the event (in general with plotFullPulse)
! plotPulseAll()
! plot pulses above threshold (in general with plotZoomPulse)
! plotPulseEvent()
! plot pulses in rows iRowT using common time for plot
! plotPulseRow(iRowT)
! plot pulses of groups iGroupV
! plotPulseVeto(iGroupV)
   PUBLIC :: plotPulseAll,plotPulseEvent,plotPulseRow,plotPulseVeto
   
   PRIVATE :: calcTMax,calcT0,check4noise
   PRIVATE :: calcSignalSum,calcSignalSlideSum,calcSignalFix
   PRIVATE :: calcSignalExp,calcSignalExpT0,amplitudeExp
   PRIVATE :: fitOnePulse,fitForT0,riseTime,fallTime
   PRIVATE :: PULSECHI,pulseShape,pulseShapePAW
   PRIVATE :: plotPulse,plotZoomPulse,plotFullPulse    
   
   REAL(KIND=minuitPrec),PRIVATE ::thisAmpl,thisT0,thisfTime,thisrTime,thisped

! default time window for pulse fitting and ploting [Tmin,Tmax] timebins
   INTEGER,PARAMETER,PRIVATE :: iTminDef =-200 & ! Tmin = T0 + iTminDef
                              , iTmaxDef = 700   ! Tmax = Tmin + iTmaxDef
   
CONTAINS

! >>>>  CALIBRATION  >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>   
   SUBROUTINE Calibrate()
      INTEGER :: iGroup  !,iRow
      REAL :: ampl,T0
      
! get the signal amplitudes
      CALL calcAmplitudeFoil()    
! fit pulses 
!      CALL fitAllPulses()  

! fill calibration histograms
      DO iGroup = 1,nGroups
         CALL getSignalGroup(iGroup,ampl,T0)
!      DO iRow=1,nRows        
!         CALL GmaxSignalRow(iRow,iGroup,ampl,T0)      
         CALL calSignal(iGroup,ampl,T0,Pedestal(iGroup))
      ENDDO
     
   END SUBROUTINE Calibrate
! <<<<  CALIBRATION  <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<

   SUBROUTINE calcPedestal()
      REAL :: ped0,ped1
      INTEGER :: nTwidth,iGroup
      PARAMETER (nTwidth=5)

! DenseData - do nothing 
      IF( nTbin==0 )RETURN
! Pedestals are read in - do nothing
      IF( .NOT.LpedestalEvent )RETURN
      
      DO iGroup = 1,nGroups
         ped0 = sum(IADC(1:nTwidth,iGroup))
         ped1 = sum(IADC(nTbin-nTwidth+1:nTbin,iGroup))
         Pedestal(iGroup) = max(ped0,ped1)/REAL(nTwidth)
      ENDDO 
   END SUBROUTINE calcPedestal
   SUBROUTINE subPedestal()
      INTEGER :: iGroup
      
! DenseData - do nothing 
      IF( nTbin==0 )RETURN

      DO iGroup = 1,nGroups
         ADC(:,iGroup) = IADC(:,iGroup) - Pedestal(iGroup)
      ENDDO
      
! Groups for which no data were read in have IADC=-1, set ADC=0
      WHERE( IADC<0 )
         ADC = 0.0
      END WHERE
   END SUBROUTINE subPedestal
   
   SUBROUTINE correctGain()
! correct the gain of each channel **********************************
      INTEGER :: iPad
      REAL :: ampl,T0
      
      DO iPad=1,nPads
         CALL GetSignal(iPad,ampl,T0)
         ampl = ampl / gainCor(PadGroup(iPad))
         CALL FillSignal(iPad,ampl,T0)
      ENDDO
   END SUBROUTINE correctGain
      
! fill the signal amplitudes for each pad, 3 routines *************
   SUBROUTINE calcAmplitudePad()
! get the signal amplitudes for each pad individually
      INTEGER :: I,iGroup,iPad,iTMax
      REAL :: amplitude,T0
      
! for dense data nTbin == 0 -> nothing to do
      IF( nTbin==0 )RETURN
            
      DO iGroup=1,nGroups
! get TMax
         CALL calcTMax(ADC(:,iGroup),iTMax)
         T0 = iTmax       
         
! some ways to calculate the signal from the pulse shape  
         CALL calcSignalFix(iGroup,iTMax,amplitude)
!         CALL calcSignalSlideSum(iGroup,-1,30,amplitude,T0)
!        CALL calcSignalSum(iGroup,iTMax,amplitude)
!         CALL calcSignalExp(iGroup,iTMax,amplitude,T0)
         
! copy this signal to all pads in the group
         DO I=1,nPadsinGroup(iGroup)
            iPad = iPadofGroup(I,iGroup)
            CALL FillSignal(iPad,amplitude,T0)
         ENDDO
      ENDDO
   END SUBROUTINE calcAmplitudePad
   
   SUBROUTINE calcAmplitudeFoilNEW()
! get the signal amplitudes setup with dispersion foil
      INTEGER :: I,iRow,iG,iGroup,iGmax,iPad,iTMax,iT0
      REAL :: amplmax,amplitude,T0,T,SL ! ,TMax
      
! for dense data nTbin == 0 -> nothing to do
      IF( nTbin==0 )RETURN
      
      DO iRow=1,nRows
         amplmax = 0.0
         iT0 = 0
         iGmax = 0
         DO iG=1,nGroupsinRow(iRow)
            iGroup = iGroupofRow(iG,iRow)

! get TMax and amplitude
            CALL calcTMax(ADC(:,iGroup),iTMax)         
            CALL calcSignalFix(iGroup,iTMax,amplitude)

            CALL check4noise(iGroup,iTMax,amplitude)
            
! get iT0 from largest signal            
            IF( amplitude > amplmax )THEN
               amplmax = amplitude
               iT0 = iTMax
               iGmax = iGroup
            ENDIF
         ENDDO

! get better T0 from quartic fit
         T0  = iT0
!         TMax = T0
         IF( iGmax /= 0 )THEN
            iT0 = max(iT0-10,1)
            CALL Qfit( ADC(iT0:iT0+24,iGmax) , T , SL)
            IF( SL > 0.0 ) T0 = T + iT0-1
         ENDIF

         DO iG=1,nGroupsinRow(iRow)
            iGroup = iGroupofRow(iG,iRow)
            
            CALL calcSignalSlideSum(iGroup,INT(T0)+35,100,amplitude,T)

!            IF( iGroup == iGmax .and. iGroup/=1 ) &
!               CALL plotZoomPulse(iGroup,amplitude,TMax)

! copy this signal to all pads in the group
            DO I=1,nPadsinGroup(iGroup)
               iPad = iPadofGroup(I,iGroup)
               CALL FillSignal(iPad,amplitude,T0)
            ENDDO
         ENDDO
      ENDDO
   END SUBROUTINE calcAmplitudeFoilNEW

   SUBROUTINE calcAmplitudeFoil()
! get the signal amplitudes setup with dispersion foil
      INTEGER :: I,iRow,iG,iGroup,iPad,iTMax,iT0,iT,nT
      REAL :: amplmax,amplpad,amplitude,T0,Tpad
      
! for dense data nTbin == 0 -> nothing to do
      IF( nTbin==0 )RETURN
      
      DO iRow=1,nRows
         amplmax = 0  
         iT0 = 0
         DO iG=1,nGroupsinRow(iRow)
            iGroup = iGroupofRow(iG,iRow)
! get TMax
            CALL calcTMax(ADC(:,iGroup),iTMax)         
! get amplitude  
            IF( iGroup<=4 )THEN
! noisy channels
               CALL calcSignalSlideSum(iGroup,-1,50,amplitude,T0)
               iTmax = T0
            ELSE
               CALL calcSignalFix(iGroup,iTMax,amplitude)
            ENDIF
! get iT0 from largest signal            
            IF( amplitude > amplmax )THEN
               amplmax = amplitude
               iT0 = iTmax
            ENDIF
         ENDDO
         
         IF( iT0<16 ) iT0=16
         
         DO iG=1,nGroupsinRow(iRow)
            iGroup = iGroupofRow(iG,iRow)
            
            amplmax = 0.0
            amplpad = 0.0
            Tpad    = 0.0
            DO iT = iT0-15,nTbin-50
               IF( iT<=iT0+15 )THEN
                  nT = 30
               ELSEIF(iT<=iT0+50) THEN
                  nT = 30 + (iT-iT0-15)*2
               ELSE
                  nT = 100
               ENDIF
               CALL calcSignalSlideSum(iGroup,iT,nT,amplitude,T0)
               IF( amplitude > amplmax )THEN
                  amplmax = amplitude
                  Tpad = T0
                  IF( iT<=iT0+50 ) amplpad = amplitude
               ENDIF
            ENDDO

! copy this signal to all pads in the group
            DO I=1,nPadsinGroup(iGroup)
               iPad = iPadofGroup(I,iGroup)
               CALL FillSignal(iPad,amplpad,TPad)
            ENDDO
         ENDDO
      ENDDO
   END SUBROUTINE calcAmplitudeFoil

   SUBROUTINE calcAmplitudeRow()
      INTEGER :: iRow,iG,iGroup,iP,iPad,iT0row
      REAL :: Amplitude,T0row,ampl,time
      
! for dense data nTbin == 0 -> nothing to do
      IF( nTbin==0 )RETURN
      
! get some independent numbers first
      CALL calcAmplitudePad()
      
! redetermine the amplitude of small pulses
! use the amplitude weighted average T0 
      RowLoop: DO iRow=1,nRows
         CALL GetSignalRow(iRow,Amplitude,T0row)
         iT0row = T0row
         
         GroupLoop: DO iG=1,nGroupsinRow(iRow)
            iGroup = iGroupofRow(iG,iRow)
            ampl = -100.5
            PadLoop: DO iP=1,nPadsinGroup(iGroup)
               iPad = iPadofGroup(iP,iGroup)
               IF( PadRow(iPad) /= iRow ) cycle PadLoop 

               CALL GetSignal(iPad,Amplitude,time)
               IF( Amplitude < smallSignal )THEN
! small pulse, redetermine amplitude
                  IF( ampl < -100.0 )THEN
! first in this group and row, calculate
!                   CALL calcSignalFix(iGroup,iT0row,ampl)
!                   CALL calcSignalSum(iGroup,iT0row,ampl)
                    CALL calcSignalExpT0(iGroup,iT0row,ampl)
                  ENDIF
! copy to pads in group
                  CALL FillSignal(iPad,ampl,T0row)
               ENDIF
            ENDDO PadLoop
         ENDDO GroupLoop
      ENDDO RowLoop
      
   END SUBROUTINE calcAmplitudeRow
   
   SUBROUTINE check4noise(iGroup,iTmax,amplitude)
! check for noisy channels
      INTEGER,INTENT(IN) :: iGroup
      INTEGER,INTENT(INOUT) :: iTmax
      REAL,INTENT(INOUT) ::amplitude
      INTEGER :: i,iT
      REAL :: ampl,aGmax

! do nothing for signals below threshold      
      IF( amplitude < minSignal ) RETURN

! if -ADC falls to less than amplitude*0.3 within 4 bins of maximum
! -> noise, recalculate amplitude with checking
      iT = iTmax
      DO i=iT-4,iT+4
         CALL calcSignalFix(iGroup,i,ampl)
         IF( ampl < amplitude*0.3 )THEN
            amplitude = 0.0
            EXIT
         ENDIF
      ENDDO
               
      IF( amplitude <= 0.0 )THEN
         aGmax = 0.0
         iTmax = 0
         DO iT=5,nTbin-5
            CALL calcSignalFix(iGroup,iT,amplitude)
            DO i=iT-4,iT+4
               CALL calcSignalFix(iGroup,i,ampl)
               IF( ampl < amplitude*0.3 )THEN
                  amplitude = 0.0
                  EXIT
               ENDIF
            ENDDO
            IF( amplitude > aGmax )THEN
               aGmax = amplitude
               iTmax = iT
            ENDIF
         ENDDO
         amplitude = aGmax
      ENDIF
   END SUBROUTINE check4noise
      
! find T0 of signal *************************************************
   SUBROUTINE calcTMax(pulse,iT0)
! time of maximum amplitude
      REAL,DIMENSION(:),INTENT(IN) :: pulse
      INTEGER,INTENT(OUT) :: iT0
!      INTEGER :: iT
!      REAL :: max
      
!      max = HUGE(max)
!      iT0 = 0
!      DO iT=1,nTbin
!         IF( pulse(iT) < max )THEN
!            max = pulse(iT)
!            iT0 = iT
!         ENDIF
!      ENDDO
 
       iT0 = minloc(pulse)
   END SUBROUTINE calcTMax
   
   SUBROUTINE calcT0(pulse,ampl,Tmax,T0)
! time of rise: pulse(T0) = ampl(Tmax)/2
      REAL,DIMENSION(:),INTENT(IN) :: pulse
      REAL,INTENT(IN) :: ampl,Tmax
      REAL,INTENT(OUT) :: T0
      INTEGER :: iTmax,iT,iT1,i,nT,ierr
      REAL :: a,b,deltaT
      REAL,DIMENSION(:),ALLOCATABLE :: x,y,ss
      
! only if there is a signal
      IF( ampl<=0.0 )THEN
         T0 = TMax
         RETURN
      ELSEIF( ampl<=smallSignal )THEN
! use default
         T0 = Tmax-Tshift 
         RETURN
      ENDIF
        
      nT = INT(Tshift)     
      allocate(x(nT), stat=ierr)
      allocate(y(nT), stat=ierr)
      allocate(ss(nT),stat=ierr)
      
       
      iTmax = min(INT(Tmax),nTbin-nT)
! get a starting point
      DO iT=iTmax-10,1,-1
         IF( abs(pulse(iT))<ampl*0.5 )EXIT
      ENDDO
      
! iT is in the middle of [iT1,iT1+nT]
      iT1 = max(0,iT-nT/2)

! linear regression around iT1: pulse(iT) = a + b*iT
      x  = (/ (real(i) , i=iT1+1,iT1+nT) /)
      y  = pulse(iT1+1:iT1+nT)
      ss = 1.0

      CALL LinearRegression(x,y,ss,a,b,ierr)    
      deallocate(x)
      deallocate(y)
      deallocate(ss)
      
! slope should be negative
      IF( ierr==0 .AND. b<0.0 )THEN      
         T0 = (-ampl*0.5 - a)/b
         deltaT = T0-Tmax
         IF( T0-Tmax < -100.0) T0 = Tmax-100.0
         IF( T0-Tmax >= 0.0 )  T0 = Tmax-Tshift
      ELSE
         T0 = Tmax-Tshift
         deltaT = 0.0
      ENDIF
      
      CALL HF_deltaT(deltaT)      
   END SUBROUTINE calcT0

! some ways to determine the signal amplitude and T0 ******************
   SUBROUTINE calcSignalFix(iGroup,iT0,ampl)
! fill the signal with T0 and the amplitude at that time
      INTEGER,INTENT(IN) :: iGroup,iT0
      REAL,INTENT(OUT) :: ampl
      
      ampl = -ADC(iT0,iGroup)
   END SUBROUTINE calcSignalFix

   SUBROUTINE calcSignalSlideSum(iGroup,iT0,nTsum,ampl,T0)
! iT0 / T0 are in the middle of the time window   
! fill the signal with T0 and the amplitude at that time
      INTEGER,INTENT(IN) :: iGroup,iT0,nTsum
      REAL,INTENT(OUT) :: ampl,T0
      INTEGER :: ifirst,ilast,nThalf
      REAL :: a,apeak

      nThalf = nTsum/2
      
      IF( iT0<=0 )THEN
! sliding sum search for largest amplitude
         a = sum(ADC(1:nTsum,iGroup))
         apeak = a
         T0 = nThalf
         do ifirst=1,nTbin-nTsum
            ilast = ifirst + nTsum
            a = a-ADC(ifirst,iGroup)+ADC(ilast,iGroup)
            IF( a<apeak )THEN
                T0 = ifirst+nThalf
                apeak = a
            ENDIF
         enddo
         ampl = -apeak/REAL(nTsum)
      ELSE
! amplitude at fixed iT0
         T0 = iT0
         ifirst = max(0,iT0-nThalf)
         ifirst = min(ifirst,nTbin-nTsum)
         ampl = -1.0*sum(ADC(ifirst+1:ifirst+nTsum,iGroup))/REAL(nTsum)
      ENDIF
   END SUBROUTINE calcSignalSlideSum
 
   SUBROUTINE calcSignalSum(iGroup,iT0,ampl)
! determine the amplitude as integral of the ADC 
! in the range iT0+iTstart to iT0+iTstart+iDeltaT
      INTEGER,INTENT(IN) :: iGroup,iT0
      REAL,INTENT(OUT) :: ampl
      INTEGER :: iTmin,iTmax,iTstart,iDeltaT
      
      iTstart = riseTime(REAL(iT0),iGroup)
      iDeltaT = iTstart + fallTime(REAL(iT0),iGroup)
! stay in range [1,nTbin]      
      iTmin = max(1,iT0-iTstart)
      iTmin = min(nTbin-iDeltaT,iTmin)
      iTmax = iTmin+iDeltaT
      ampl = -1.0 * sum( ADC(iTmin:iTmax,iGroup) )/REAL(iDeltaT+1)
      
   END SUBROUTINE calcSignalSum
 
   SUBROUTINE calcSignalExp(iGroup,iT0,ampl,T0)
      INTEGER,INTENT(IN) :: iGroup,iT0
      REAL,INTENT(OUT) :: ampl,T0
      REAL :: TMax
      
      CALL amplitudeExp(iGroup,iT0,ampl)
      TMax = REAL(iT0)
      
! find T0 where ADC(T0) = ampl/2
      IF( ampl>0.0 ) THEN
         CALL calcT0(ADC(:,iGroup),ampl,TMax,T0) 
         CALL HF_deltaT(T0-TMax)
         ampl = ampl * exp((TMax-T0)/fallTime(T0,iGroup))
      ELSE
         T0 = 0.0
      ENDIF
   END SUBROUTINE calcSignalExp
   
   SUBROUTINE calcSignalExpT0(iGroup,iT0,ampl)
      INTEGER,INTENT(IN) :: iGroup,iT0
      REAL,INTENT(OUT) :: ampl
      INTEGER :: iTmax
      REAL :: T0
         
      iTmax = iT0 + Tshift
      CALL amplitudeExp(iGroup,iTmax,ampl)
      T0 = REAL(T0)
      
! correct from TMax to T0
      ampl = ampl * exp(REAL(Tshift)/fallTime(T0,iGroup))

   END SUBROUTINE calcSignalExpT0

   SUBROUTINE amplitudeExp(iGroup,iT0,Amplitude)
! calculates the amplitude at iT0 from the average 
! in the range iT0+iTstart to iT0+iTstart+iDeltaT
! assuming an exponential with fallTime 
      INTEGER,INTENT(IN) :: iGroup,iT0
      REAL,INTENT(OUT) :: Amplitude
      INTEGER :: iG,iT,iDT,IERR
! iDeltaTmin should NEVER be larger than iDeltaT
      INTEGER,PARAMETER :: iTstart=50,iDeltaT=300,iDeltaTmin=100
      REAL :: ampl,scaleF
      REAL,ALLOCATABLE,SAVE,DIMENSION(:) :: scaleFac 
      LOGICAL,SAVE :: first=.TRUE.
      
! at first call calculate scaleFac      
      IF( first )THEN 
         ALLOCATE ( scaleFac(nGroups),STAT=IERR )
         scaleFac = 0.0
         DO iG = 1,nGroups
         DO iT=1,iDeltaT
            scaleFac(iG) = scaleFac(iG) &
                         + exp(-REAL(iT+iTstart)/fallTime(REAL(nTbin)/2.0,iG))
! in case of time dependent fall time use typical time: nTbin/2
         ENDDO
         ENDDO
         first = .FALSE.
      ENDIF

! plain average is equivalent to error weighted mean with error = 1ADC count
      IF( iT0+iTstart+iDeltaT <= nTbin )THEN
         ampl = sum( ADC(iT0+iTstart:iT0+iTstart+iDeltaT,iGroup) )
         ampl = -1.0 * ampl/scaleFac(iGroup) 
      ELSEIF(iT0+iTstart+iDeltaTmin <= nTbin )THEN
! recalculate scaleFac
         iDT = nTbin-iT0-iTstart
         scaleF = 0.0
         DO iT=1,iDT
            scaleF = scaleF &
                   + exp(-REAL(iT+iTstart)/fallTime(REAL(iTstart),iGroup))
         ENDDO
         ampl = sum( ADC(iT0+iTstart:iT0+iTstart+iDT,iGroup) )
         ampl = -1.0 * ampl/scaleF 
      ELSE
! pulse too late
         ampl = 0.0
      ENDIF

      Amplitude = ampl

   END SUBROUTINE amplitudeExp
      
   SUBROUTINE fitAllPulses()
      INTEGER :: iGroup,iP,iPad,iPmax
      REAL :: minAmpl,amplmax,ampl,T0

      IF( Lcalib )THEN
! for calibration only big pulses
         minAmpl = minSignalCalib
      ELSE
         minAmpl = smallSignal
      ENDIF
      
      DO iGroup = 1,nGroups
! fit the pad with largest amplitude
         iPmax = 0
         amplmax = 0.0
         DO iP = 1,nPadsinGroup(iGroup)
            iPad = iPadofGroup(iP,iGroup)
            CALL GetSignal(iPad,ampl,T0)
            IF( ampl>amplmax )THEN
               amplmax = ampl
               iPmax = iPad
            ENDIF
         ENDDO 
         
         IF( iPmax == 0 )CYCLE
         
         iPad = iPmax
         CALL GetSignal(iPad,ampl,T0)
         IF(ampl > minAmpl )THEN
            CALL fitOnePulse(iPad,ampl,T0)
            
! copy the signal to other pads in group
            DO iP = 1,nPadsinGroup(iGroup)
               iPad = iPadofGroup(iP,iGroup)
               CALL FillSignal(iPad,ampl,T0)
            ENDDO 
         ENDIF
      ENDDO
   END SUBROUTINE fitAllPulses

   SUBROUTINE fitOnePulse(iPad,ampl,T0)
      INTEGER , INTENT(IN) :: iPad
      REAL , INTENT(INOUT) :: ampl,T0
      
      CHARACTER(len=80) :: CHTITL
      INTEGER,PARAMETER :: NFITP=9
      INTEGER :: IPRT,IERR,iGroup
      REAL :: rTime,fTime,ped,amplfit,T0fit
      REAL,DIMENSION(100) :: Xdata
      CHARACTER(len=10),DIMENSION(NFITP) :: PCH
      REAL(KIND=minuitPrec),DIMENSION(NFITP) :: PVAL,PERR,BND1,BND2
      REAL(KIND=minuitPrec),DIMENSION(NFITP,NFITP) :: PMAT
      REAL(KIND=minuitPrec) :: CHI2
      
      iGroup  = PadGroup(iPad)
      rTime   = riseTime(T0,iGroup)
      fTime   = fallTime(T0,iGroup)
      ped     = 0.0 ! fit pedestal subtracted ADC

! names
      PCH(1)  = "Amplitude "
      PCH(2)  = "T0        "
      PCH(3)  = "Risetime  "
      PCH(4)  = "Falltime  "
      PCH(5)  = "Pedestal  "
      PCH(6)  = "T0Start   "
      PCH(7)  = "iGroup    "
      PCH(8)  = "iTmin     "
      PCH(9)  = "iTmax     "
! start values
      PVAL(1) = ampl
      PVAL(2) = T0
      PVAL(3) = rTime
      PVAL(4) = fTime
      PVAL(5) = ped
      PVAL(6) = T0
      PVAL(7) = iGroup
      PVAL(8) = 0.0
      PVAL(9) = 0.0
      
! error estimates
      PERR(1) = ampl*0.1+5.0
      PERR(2) = 5.0
      IF( LCALIB )THEN
         PERR(3) = 2.0
         PERR(4) = 20.0
         PERR(5) = 3.0
      ELSE
         PERR(3) = -1.0
         PERR(4) = -1.0
         PERR(5) = -1.0
      ENDIF
      PERR(6) = -1.0
      PERR(7) = -1.0
      PERR(8) = -1.0
      PERR(9) = -1.0

! normaly no bounds
      BND1 = 0.0      
      BND2 = 0.0  
! MINUIT needs some help
      BND1(2) = MAX(0.0,T0-200.0)    
      BND2(2) = BND1(2)+400.0    
      IF( LCALIB )THEN
         BND1(3) = 3.5
         BND2(3) = 20.0
         BND1(4) = 200.0
         BND2(4) = 600.0
         BND1(5) = -30.0    
         BND2(5) =  30.0  
      ENDIF
              
      WRITE(UNIT=CHTITL,FMT=*)"+++ Pulse fit",iGroup
      IPRT = -1
      CALL MINUIT(PULSECHI,CHTITL,NFITP,PCH,PVAL,PERR,PMAT, &
                  BND1,BND2,CHI2,IPRT,IERR)
      
! plot pulse with start values
      ped = Pedestal(iGroup)
      IF( .NOT.Lcalib ) CALL plotPulse(iGroup,ampl,T0,rTime,fTime,ped,0.0)

! get fit results  
! PVAL(1) = amplitude at T0start -> need corection       
      amplfit  = abs(PVAL(1))*exp((PVAL(6)-PVAL(2))/abs(PVAL(4)))
      T0fit    = PVAL(2)
      rTime = abs(PVAL(3))
      fTime = abs(PVAL(4))
      ped   = PVAL(5)+Pedestal(iGroup)
      CALL plotPulse(iGroup,amplfit,T0fit,rTime,fTime,ped,REAL(CHI2))

      IF( .NOT.Lcalib )THEN
! ntuple comparing start and fit values
         Xdata(1) = eventNumber
         Xdata(2) = iGroup
         Xdata(3) = ampl
         Xdata(4) = amplfit
         Xdata(5) = T0
         Xdata(6) = T0fit
         Xdata(7) = CHI2
         CALL NF_pulsefit(Xdata)
      ENDIF


! don't change the amplitude and T0 if fit failed
! failed: VERY bad CHI2
! OK: VERY good CHI2, reasonable IERR
      IF( CHI2<10000.0 .AND. &
        (CHI2<700.0 .OR. IERR == 0 .OR. IERR == 1 .OR. IERR == 2 ) )THEN

         ampl = amplfit
         T0 = T0fit

! fill the pedestal, risetime and falltime in histograms     
         IF( Lcalib ) CALL calRFtime(iGroup,rTime,fTime,T0,ped)
      ENDIF
      
   END SUBROUTINE fitOnePulse
   SUBROUTINE fitForT0(iGroup,ampl,TMax,T0)
      INTEGER , INTENT(IN) :: iGroup
      REAL , INTENT(IN) :: ampl,TMax
      REAL , INTENT(OUT) :: T0
      
      CHARACTER(len=80) :: CHTITL
      INTEGER,PARAMETER :: NFITP=9
      INTEGER :: IPRT,IERR
      CHARACTER(len=10),DIMENSION(NFITP) :: PCH
      REAL(KIND=minuitPrec),DIMENSION(NFITP) :: PVAL,PERR,BND1,BND2
      REAL(KIND=minuitPrec),DIMENSION(NFITP,NFITP) :: PMAT
      REAL(KIND=minuitPrec) :: CHI2
      
! only if there is a signal
      IF( ampl<=0.0 )THEN
         T0 = TMax
         RETURN
      ENDIF
      
! names
      PCH(1)  = "Amplitude "
      PCH(2)  = "T0        "
      PCH(3)  = "Risetime  "
      PCH(4)  = "Falltime  "
      PCH(5)  = "Pedestal  "
      PCH(6)  = "T0Start   "
      PCH(7)  = "iGroup    "
      PCH(8)  = "iTmin     "
      PCH(9)  = "iTmax     "
! start values
      PVAL(1) = ampl
      PVAL(2) = TMax-Tshift
      PVAL(3) = riseTime(TMax,iGroup)
      PVAL(4) = fallTime(TMax,iGroup)
      PVAL(5) = Pedestal(iGroup)
      PVAL(6) = TMax
      PVAL(7) = iGroup
      PVAL(8) = max(1.0,TMax-100.0)
      PVAL(9) = TMax+10.0
! error estimates
      PERR(1) = -1.0
      PERR(2) =  5.0
      PERR(3) = -1.0
      PERR(4) = -1.0
      PERR(5) = -1.0
      PERR(6) = -1.0
      PERR(7) = -1.0
      PERR(8) = -1.0
      PERR(9) = -1.0

! normaly no bounds
      BND1 = 0.0      
      BND2 = 0.0  
! MINUIT needs some help
      BND1(2) = TMax-100.0    
      BND2(2) = TMax    
              
      WRITE(UNIT=CHTITL,FMT=*)"+++ T0 fit",iGroup
      IPRT = -1
      CALL MINUIT(PULSECHI,CHTITL,NFITP,PCH,PVAL,PERR,PMAT, &
                  BND1,BND2,CHI2,IPRT,IERR)
      
      T0 = PVAL(2)
      
   END SUBROUTINE fitForT0

   FUNCTION riseTime(T0,iGroup) RESULT(rTime)
! returns rise time depending on provided information   
      INTEGER,INTENT(IN) :: iGroup
      REAL,INTENT(IN) :: T0
      REAL :: rTime
      
      IF( size(riseTimeGr)>0 )THEN
! TPCrtime.txt : depending on readout channel
         rTime = riseTimeGr(iGroup)
      ELSEIF( size(riseTimePar)>0 .AND. T0>0.0 )THEN
! Threshold.txt, polynomial : depending on time bin
         rTime = POLY(T0,riseTimePar)
      ELSE
! default
         rTime = Tshift
      ENDIF            
   END FUNCTION riseTime

   FUNCTION fallTime(T0,iGroup) RESULT(fTime)
! returns fall time depending on provided information   
      INTEGER,INTENT(IN) :: iGroup
      REAL,INTENT(IN) :: T0
      REAL :: fTime
      
      IF( size(fallTimeGr)>0 )THEN
! TPCftime.txt : depending on readout channel
         fTime = fallTimeGr(iGroup)
      ELSEIF( size(fallTimePar)>0 .AND. T0>0.0 )THEN
! Threshold.txt, polynomial : depending on time bin
         fTime = POLY(T0,fallTimePar)
      ELSE
! default
         fTime = Tshift
      ENDIF            
   END FUNCTION fallTime
   
   SUBROUTINE PULSECHI(NPAR,GRAD,FVAL,XVAL,IFLAG)
! subroutine to be called by MINUIT - calculates the chisq
      INTEGER,INTENT(IN) :: IFLAG , NPAR
      INTEGER,PARAMETER :: NFITP=9
      REAL(KIND=minuitPrec),DIMENSION(NFITP),INTENT(IN) :: XVAL
      REAL(KIND=minuitPrec),DIMENSION(NFITP),INTENT(OUT) :: GRAD
      REAL(KIND=minuitPrec),INTENT(OUT) :: FVAL
      REAL(KIND=minuitPrec) :: CHISQ,diff
      REAL(KIND=minuitPrec) :: T,ampl,T0,rTime,fTime,ped,T0start
      INTEGER :: iGroup,iT,iTmin,iTmax,myPar
   
! dummy to avoid compiler warnings
      IF( IFLAG==9 ) myPar=NPAR

      CHISQ = 0.0_minuitPrec
      GRAD  = 0.0_minuitPrec
      T0    = XVAL(2)
      rTime = ABS(XVAL(3))
      fTime = ABS(XVAL(4))
      ped   = XVAL(5)
      T0start = XVAL(6)
      iGroup  = anint(XVAL(7))
      iTmin  = anint(XVAL(8))
      iTmax  = anint(XVAL(9))
! this is the amplitude at fixed T0start
! -> XVAL(1) and XVAL(2) should be uncorrelated
      ampl  = ABS(XVAL(1))/exp((T0start-T0)/fTime)
      
      IF( iTmax == 0 )THEN
! window for fit: use defaults iTmaxDef bins, iTminDef bins wrt T0, within [1,nTbin]
         iTmin = T0 + iTminDef
         IF( iTmin<1 ) iTmin=1
         IF( iTmin>nTbin-iTmaxDef ) iTmin=nTbin-iTmaxDef
         iTmax = iTmin + iTmaxDef
      ENDIF

      DO iT = iTmin,iTmax
         T = iT
! underflow, don't include in fit (or IADC not filled)
         IF( IADC(iT,iGroup) <= 1 )cycle
         diff  = (ADC(iT,iGroup)-pulseShape(T,ampl,T0,rTime,fTime,ped)) 
         CHISQ = CHISQ + diff*diff
      ENDDO
      
      FVAL = CHISQ
   END SUBROUTINE PULSECHI
   FUNCTION pulseShape(T,ampl,T0,rTime,fTime,ped) RESULT(A)
! 'fast' version of pulseShapePAW
! apply any change to both pulseShape and pulseShapePAW
      REAL(KIND=minuitPrec),INTENT(IN) :: T,ampl,T0,rTime,fTime,ped
      REAL(KIND=minuitPrec) :: A,x
      
      x = T-T0
      
      A = ped - ampl / &
          (1.0_minuitPrec+exp(-x/rTime)) / exp(x/fTime)
      
   END FUNCTION pulseShape
   
   FUNCTION pulseShapePAW(T) RESULT(A)
! keep T and A REAL so function can be used in call to HFUNC
! other parameters are passed as global variables
      REAL,INTENT(IN) :: T
      REAL :: A
      REAL(KIND=minuitPrec) :: x,rTime,fTime,xrTime,xfTime,ampl
      
! this is an empirical offset 
      x = T-thisT0
      ampl = thisAmpl
      
      rTime = thisrTime
      fTime = thisfTime
! try to avoid arithmetic exceptions
      IF( abs(rTime)<1.0E-8_minuitPrec ) rTime=1.0E-8_minuitPrec
      IF( abs(fTime)<1.0E-8_minuitPrec ) fTime=1.0E-8_minuitPrec
      xrTime = x/rTime
      xfTime = x/fTime
      IF( xrTime >  20.0_minuitPrec) xrTime = 20.0_minuitPrec
      IF( xfTime >  40.0_minuitPrec) xfTime = 40.0_minuitPrec
      IF( xrTime < -40.0_minuitPrec)THEN
          xrTime = -40.0_minuitPrec
          xfTime = xrTime*rTime/fTime
      ENDIF
      
      A = thisped - &
          ampl / (1.0_minuitPrec+exp(-xrTime)) / exp(xfTime)
      
   END FUNCTION pulseShapePAW
   
   SUBROUTINE plotPulse(iGroup,ampl,T0,rTime,fTime,ped,chisq)
      INTEGER,INTENT(IN) :: iGroup
      REAL,INTENT(IN) :: ampl,T0,fTime,rTime,ped,chisq
      INTEGER :: iT0,it,i,iHist,nbin,ioff,nfac
      REAL :: x,pmin,pmax
      INTEGER,PARAMETER :: nbmax=200
      REAL,DIMENSION(nbmax) :: pulse
      CHARACTER(len=80) :: CHtitle
      
      IF( nTbin==0 )RETURN
      
! set time window for plot (same as in PULSECHI),
! average over nfac time bins
      nfac = 5
      nbin = iTmaxDef-iTminDef
      IF( nbin/nbmax > nfac )  nfac = nbin/nbmax
      nbin = nbin/nfac         
      ioff = iTminDef/nfac
      
      iT0 = anint(t0) 
      DO i=1,nbin
! T0 is at i=ioff, nfac bins averaged
         it = (i+ioff-1)*nfac + iT0 
         x = 0.0
         IF( it>=0 .AND. it<=nTbin-nfac )THEN
            x = sum( IADC(it+1:it+nfac,iGroup) )
            pulse(i) = x/REAL(nfac)
         ELSE
            pulse(i) = 0.0
         ENDIF
      ENDDO
      
! fill pulse in histogram
      WRITE(UNIT=CHtitle,FMT="(A,I7,A,I4,A,F8.2,A,F8.2)") &
          "Event",EventNumber," Group",iGroup , &
          "    amplitude",ampl,"   chisq",chisq          
      pmin = T0+ioff*nfac
      pmax = T0+(nbin+ioff)*nfac
      CALL HF_pulse(CHtitle,pulse(1:nbin),pmin,pmax,iHist)
      
! add function
      IF( ampl > 0.0 .AND. iHist > 0 )THEN
         thisAmpl = ampl
         thisT0   = T0
         thisrTime = rTime 
         thisfTime = fTime
         thisped  = ped
         CALL HFUNC(iHist,pulseShapePAW)
      ENDIF
   END SUBROUTINE plotPulse
   SUBROUTINE plotFullPulse(iGroup,ampl,T0)
      INTEGER,INTENT(IN) :: iGroup
      REAL,INTENT(IN) :: ampl,T0
      INTEGER :: it,i,nfac,nTb
      INTEGER,SAVE :: iHist
      REAL :: x
      INTEGER,PARAMETER :: nbin=200
      REAL,DIMENSION(nbin) :: pulse
      CHARACTER(len=80) :: CHtitle
      
      nTb = nTbGroup(iGroup)
      IF( nTb<=1 )RETURN
      
      IF( nTb > nbin)THEN
         nfac = nTb/nbin
!        nfac = 5
         nTb = nbin
      ELSE
         nfac = 1
      ENDIF
      
      DO i=1,nTb
         it = (i-1)*nfac  
         x = sum( IADC(it+1:it+nfac,iGroup) )
         pulse(i) = x/REAL(nfac)
      ENDDO
      
! fill pulse in histogram
      WRITE(UNIT=CHtitle,FMT="(A,I7,A,I4,A,F8.2,A,F8.1)") &
          "Event",EventNumber," Group",iGroup &
         ,"    amplitude",ampl,"   T0",T0
!      CALL HF_pulse(CHtitle,pulse,0.0,REAL(nTbin),iHist)
      CALL HF_pulse(CHtitle,pulse(1:nTb),T0,T0+REAL(nTb*nfac),iHist)
!      CALL HF_pulse(CHtitle,pulse(1:nTb),0.0,REAL(nTb*nfac),iHist)
!      CALL HF_pulse(CHtitle,ADC(1:100,iGroup),0.0,REAL(nbin),iHist)
      
   END SUBROUTINE plotFullPulse
   SUBROUTINE plotZoomPulse(iGroup,ampl,T0)
      INTEGER,INTENT(IN) :: iGroup
      REAL,INTENT(IN) :: ampl,T0
      INTEGER :: it,i,nfac,iHist,ifirst
      REAL :: x
!      INTEGER,PARAMETER :: nbin=200,ntot=200,noffset=50
      INTEGER,PARAMETER :: nbin=200,ntot=1000,noffset=200
      REAL,DIMENSION(nbin) :: pulse
      CHARACTER(len=80) :: CHtitle
      
      IF( nTbin==0 )RETURN
      
      nfac = ntot/nbin
      
      ifirst = min( int(T0)-noffset,nTbin-ntot ) -1
      if( ifirst<0 ) ifirst=0
      
      DO i=1,nbin
         it = (i-1)*nfac + ifirst
         x = sum( ADC(it+1:it+nfac,iGroup) )
         pulse(i) = x/REAL(nfac)
      ENDDO
      
! fill pulse in histogram
      WRITE(UNIT=CHtitle,FMT="(A,I7,A,I4,A,F8.2,A,F8.1)") &
          "Event",EventNumber," Group",iGroup &
         ,"    amplitude",ampl,"   T0",T0
      CALL HF_pulse(CHtitle,pulse,REAL(ifirst),REAL(ifirst+ntot),iHist)
      
   END SUBROUTINE plotZoomPulse
   SUBROUTINE plotPulseAll()
      INTEGER :: iGroup
      REAL :: ampl, T0
      
      DO iGroup = 1,nGroups
         CALL GetSignalGroup(iGroup,ampl,T0)
         
         CALL plotFullPulse(iGroup,ampl,T0)
!         CALL plotZoomPulse(iGroup,ampl,T0)
      ENDDO
   END SUBROUTINE plotPulseAll
   SUBROUTINE plotPulseEvent()
      INTEGER :: iGroup
      REAL :: ampl, T0
      
      DO iGroup = 1,nGroups
! need proper time for window
         CALL GetSignalGroup(iGroup,ampl,T0)
         IF( ampl>0.0 ) &
         CALL plotZoomPulse(iGroup,ampl,T0)
      ENDDO
   END SUBROUTINE plotPulseEvent
   SUBROUTINE plotPulseRow(iRowT)
      INTEGER,DIMENSION(:),INTENT(IN) :: iRowT
      INTEGER :: iR,iRow,iG,iGroup,iGmax
      REAL :: ampl, T0, Trow
      
      DO iR = 1,size(iRowT)
         iRow = iRowT(iR)
! need proper time for window
         CALL GmaxSignalRow(iRow,iGmax,ampl,Trow)
         DO iG = 1,nGroupsinRow(iRow)
            iGroup = iGroupofRow(iG,iRow) 
            CALL GetSignalGroup(iGroup,ampl,T0)
            IF( ampl>0.0 ) &
            CALL plotZoomPulse(iGroup,ampl,Trow)
         ENDDO
      ENDDO
   END SUBROUTINE plotPulseRow
   SUBROUTINE plotPulseVeto(iGroupV)
      INTEGER,DIMENSION(:),INTENT(IN) :: iGroupV
      INTEGER :: iG, iGroup
      REAL :: ampl, T0
      
      DO iG = 1,size(iGroupV)
         iGroup = iGroupV(iG)
! need proper time for window
         CALL GetSignalGroup(iGroup,ampl,T0)
         CALL plotFullPulse(iGroup,ampl,T0)
      ENDDO
   END SUBROUTINE plotPulseVeto
END MODULE TPCcalcAmplitude