MODULE LCIO4TPC
#define PTRTYPE integer
   USE LCIOapi
   USE TPCdata
   PRIVATE
   
! OpenReadLCIO(CHFILE,IOS)
! ReadLCIOnTbin(nTbin)
! ReadLCIOComment(chComment,chDate)
! CloseReadLCIO()
   PUBLIC :: OpenReadLCIO,CloseReadLCIO,ReadLCIOComment,ReadLCIOnTbin

! ReadLCIONextEvent(IERR)
   PUBLIC :: ReadLCIONextEvent

! OpenWriteLCIO()
! CloseWriteLCIO()
! WriteLCIOHeader(chComment,chDate)
! WriteLCIOEvent(iflag)
   PUBLIC :: OpenWriteLCIO,CloseWriteLCIO,WriteLCIOHeader,WriteLCIOEvent

   PRIVATE :: ReadLCIOTPChit,ReadLCIOMCtrack

   PTRTYPE,PRIVATE :: reader, writer

CONTAINS
   SUBROUTINE OpenReadLCIO(CHFILE,IOS)
      CHARACTER(LEN=*),INTENT(IN) :: CHFILE
      INTEGER,INTENT(OUT) :: IOS
      INTEGER :: status
      LOGICAL,SAVE :: first=.TRUE.
      
      IF( first )THEN
         reader = lcrdrcreate()
         first = .FALSE.
      ELSE
         status = lcrdrclose( reader )
      ENDIF
      
      status = lcrdropen( reader, CHFILE )
      IF( status == LCIO_SUCCESS )THEN
         IOS = 0
      ELSE
         IOS = 1
      ENDIF
   END SUBROUTINE OpenReadLCIO
   SUBROUTINE ReadLCIOnTbin(nTbin)
      INTEGER,INTENT(OUT) :: nTbin
      PTRTYPE :: event,tphcol,tphit
      INTEGER :: nTb,iHit,nHits,ievent,iflag
      LOGICAL,SAVE :: first=.TRUE.
      
      IF( .not.first )THEN
         nTbin = -1
         RETURN
      ENDIF 
      
      nTbin = -2
      first = .FALSE.
      
      event_loop: DO ievent=1,10
         event  = lcrdrreadnextevent(reader,LCIO_READ_ONLY) 
         IF( event==0 ) cycle event_loop
         
         tphcol = lcevtgetcollection(event,LCIO_TPCHIT) 
         IF( tphcol==0 ) cycle event_loop
         
         iflag = lccolgetflag(tphcol)
         IF( .NOT.BTEST(iflag,LCIO_TPCBIT_RAW) )THEN
! no raw data stored
            nTbin = 0
            RETURN
         ENDIF
         
         nHits = lccolgetnumberofelements( tphcol )
       
         hit_loop: DO iHit=1,nHits
            tphit = lccolgetelementat(tphcol,iHit)
            IF(tphit==0) cycle hit_loop
            nTb  = lctphgetnrawdatawords(tphit)
            IF( nTb > 0 )THEN
               nTbin = nTb
               exit event_loop
            ENDIF
         ENDDO hit_loop
         
      ENDDO event_loop
      
   END SUBROUTINE ReadLCIOnTbin
   SUBROUTINE ReadLCIOComment(chComment,chDate)
      CHARACTER(LEN=*),INTENT(OUT) :: chComment,chDate   
      PTRTYPE :: runhdr
      
      runhdr = lcrdrreadnextrunheader( reader, LCIO_READ_ONLY )
         
      chComment = lcrhdgetdetectorname( runHdr )
      chDate = lcrhdgetdescription( runHdr )
      
   END SUBROUTINE ReadLCIOComment

   SUBROUTINE ReadLCIONextEvent(IERR)   
      INTEGER,INTENT(OUT) :: IERR
      PTRTYPE :: event
      
      IERR = 0
      
      event = lcrdrreadnextevent( reader,LCIO_READ_ONLY )
      IF( event==0 ) THEN
         IERR = -1
         RETURN
      ENDIF
      
      RunNumber   = lcevtgetrunnumber(event) 
      EventNumber = lcevtgeteventnumber(event) 
   
      CALL ReadLCIOTPChit(event)
      
      CALL ReadLCIOMCtrack(event)
      
   END SUBROUTINE ReadLCIONextEvent   
 
   SUBROUTINE ReadLCIOTPChit(event)   
      PTRTYPE,INTENT(IN) :: event
      PTRTYPE :: tphcol,tphit
      INTEGER :: iflag,nHits,iHit,iGroup,iP,iPad,nTb,i,rawData
      REAL :: Amplitude,Time0
      
      tphcol = lcevtgetcollection(event, LCIO_TPCHIT )
      IF( tphcol==0 ) RETURN
       
      nHits = lccolgetnumberofelements( tphcol )
      
      iflag = lccolgetflag(tphcol)
      IF(  nTbin>0 .AND. BTEST(iflag,LCIO_TPCBIT_RAW)&
                   .AND. .NOT.BTEST(iflag,1) )THEN
! nTbin ==0: no space for RawData
! BTEST(iflag,LCIO_TPCBIT_RAW): RawData not in LCIO
! BTEST(iflag,1): use Amplitude / Time0 instead
         
         IADC = -1

         DO iHit=1,nHits
            tphit = lccolgetelementat(tphcol,iHit)
            IF( tphit==0 ) cycle
             
            iGroup = lctphgetcellid(tphit)                   
            nTb  = lctphgetnrawdatawords(tphit)
            IF( nTb>0 .AND. nTb/=nTbin )THEN
               PRINT*,"WARNING, ReadLCIOEvent> time bins of this channel doesnt match nTbin"
               PRINT*,"EventNumber:",EventNumber," Group:",iGroup," nTb:",nTb," nTbin:",nTbin
               IF( nTb>nTbin ) nTb=nTbin
            ENDIF
            DO i=1,nTb
! LCIO feature: timebins run from [0,nTb-1]
               rawData = lctphgetrawdataword(tphit,i-1)
               IF( rawData < 0 )THEN
                  PRINT*,"WARNING, ReadLCIOEvent> IADC is negative"
                  rawData = 0
               ELSEIF( rawData > huge(IADC(1,1)) )THEN
                  PRINT*,"WARNING, ReadLCIOEvent> IADC >",huge(IADC(1,1))
                  PRINT*,"What FADC do you have? Increase kind for IADC."
                  rawData = huge(IADC(1,1))
               ENDIF
               IADC(i,iGroup) = rawData
            ENDDO
! in case nTb<nTbin
            DO i=nTb+1,nTbin
               IADC(i,iGroup) = -1
            ENDDO
         ENDDO
      ELSE
! Amplitude / Time0         
         DO iPad=1,nPads
            CALL FillSignal(iPad,0.0,1.0)
         ENDDO
         
         DO iHit=1,nHits
            tphit = lccolgetelementat(tphcol,iHit)
            IF( tphit==0 ) cycle
             
            iGroup = lctphgetcellid(tphit)
            Amplitude = lctphgcharge(tphit)
            Time0 = lctphgettime(tphit)
            DO iP=1,nPadsinGroup(iGroup)
               iPad = iPadofGroup(iP,iGroup)
               CALL FillSignal(iPad,Amplitude,Time0)
            ENDDO
         ENDDO     
      ENDIF
      
   END SUBROUTINE ReadLCIOTPChit
   
   SUBROUTINE ReadLCIOMCtrack(event)
      PTRTYPE,INTENT(IN) :: event
      PTRTYPE :: mccol,ptrack
      INTEGER :: nTrack,i

      mccol = lcevtgetcollection(event, "MCtrack" )
      IF( mccol==0 ) RETURN
      
      nTrack = lccolgetnumberofelements( mccol )
      IF( nTrack==0 ) RETURN
      
      ptrack = lccolgetelementat(mccol,1)
      IF( ptrack==0 ) RETURN
      
      IF( lcfvcgetlength(ptrack) /= 4 )THEN
         PRINT*,"Strange vector for MCtrack"
         PRINT*,(lcsvcgetfloatat(ptrack,i),i=1,lcfvcgetlength(ptrack))
         RETURN
      ENDIF
      
      MCx0    = lcsvcgetfloatat(ptrack,1)
      MCphi   = lcsvcgetfloatat(ptrack,2)
      MCz0    = lcsvcgetfloatat(ptrack,3)
      MCtheta = lcsvcgetfloatat(ptrack,4) 
       
   END SUBROUTINE ReadLCIOMCtrack
     
   SUBROUTINE CloseReadLCIO()   
      INTEGER :: status
      
      status = lcrdrclose( reader )
      status = lcrdrdelete( reader )
   END SUBROUTINE CloseReadLCIO

! -------------------------------------------------------------------
   
   SUBROUTINE OpenWriteLCIO() 
      INTEGER :: status
      
      writer = lcwrtcreate()
      
      status = lcwrtopen( writer, "thisrun.slcio",LCIO_WRITE_NEW )
      if( status==LCIO_SUCCESS) then
         print*, " opened file .slcio "
      else
         Print*, " could not open file .slcio "
      endif
   END SUBROUTINE OpenWriteLCIO
   SUBROUTINE CloseWriteLCIO() 
      INTEGER :: status
      
      status = lcwrtclose( writer )
      status = lcwrtdelete(writer)
   END SUBROUTINE CloseWriteLCIO
   
   SUBROUTINE WriteLCIOHeader(chComment,chDate)
      CHARACTER(LEN=*),INTENT(IN) :: chComment,chDate   
      PTRTYPE :: runhdr
      INTEGER :: status

      runhdr = lcrhdcreate()
      status = lcrhdsetrunnumber( runhdr, RunNumber )
      status = lcrhdsetdetectorname( runhdr, chComment  )
      status = lcrhdsetdescription( runhdr,chDate)
      status = lcwrtwriterunheader( writer, runhdr )
      status = lcrhddelete(runhdr)
   END SUBROUTINE WriteLCIOHeader
   SUBROUTINE WriteLCIOEvent(iflag)
! iflag describes whether to write RawData and/or DenseData
!     iflag = IBSET(iflag,0) ! RawData   
!     iflag = IBSET(iflag,1) ! DenseData 
      INTEGER, INTENT(IN) :: iflag
      PTRTYPE :: event,tphcol,tphit
      INTEGER :: lcioflag,status,iGroup
      REAL :: Amplitude,T0
      
      event  = lcevtcreate()
      tphcol = lccolcreate( LCIO_TPCHIT )

      status = lcevtsetrunnumber( event, RunNumber ) 
      status = lcevtseteventnumber( event,  EventNumber ) 
      
      lcioflag = 0
! write RawData
      IF( BTEST(iflag,0) ) lcioflag = IBSET(lcioflag,LCIO_TPCBIT_RAW)
! write Amplitude / Time0
      IF( BTEST(iflag,1) ) lcioflag = IBSET(lcioflag,1)

      status = lccolsetflag( tphcol, lcioflag )

      DO iGroup=1,nGroups
         tphit  = lctphcreate() 
         status = lctphsetcellid( tphit, iGroup )          

         IF( BTEST(iflag,0) )THEN
            status = lctphsetrawdata( tphit, INT(IADC(:,iGroup)), nTbin ) 
         ENDIF
         IF( BTEST(iflag,1) )THEN
            CALL getSignalGroup(iGroup,Amplitude,T0)
            status = lctphsetcharge(tphit,Amplitude)
            status = lctphsettime(tphit,T0)
         ENDIF
         status = lccoladdelement(tphcol,tphit) 
      ENDDO
 
!---- now add everything to the event
      status = lcevtaddcollection(event,tphcol,LCIO_TPCHIT)  
      status = lcwrtwriteevent(writer,event)
            
      IF( status == LCIO_ERROR )THEN
         PRINT*, "TPCTEST: error - can not write event", EventNumber
         STOP
      ENDIF
       
      status = lcevtdelete(event)
   END SUBROUTINE WriteLCIOEvent

END MODULE LCIO4TPC