*
* $Id: gputrd.F,v 3.1.1.1 2002/07/11 16:02:01 cafagna Exp $
*
* $Log: gputrd.F,v $
* Revision 3.1.1.1  2002/07/11 16:02:01  cafagna
* First GPAMELA release on CVS
*
*
*CMZ :  3.00/00 07/02/2002  12.41.53  by  Unknown
*CMZ :  2.01/00 05/04/2000  14.35.18  by  Marialuigia Ambriola
*CMZ :  2.00/00 29/02/2000  10.23.24  by  Francesco Cafagna
*CMZ :  1.02/00 09/02/2000  17.22.16  by  Francesco Cafagna
*CMZ :  1.00/02 04/04/96  16.47.16  by  Francesco Cafagna
*-- Author :    Francesco Cafagna   21/02/96
      SUBROUTINE GPUTRD(IACT,IS,ID,IPA,IT,DELOSS,PATH,
     +                  TRAPAR,NUMVOL,IMEC)
************************************************************************
*                                                                      *
* To collect and store the hit structure for TRD                       *
*                                                                      *
* Variables definition:                                                *
* IN:                                                                  *
*  IACT, integer specifing the action to be taken. It is the INWVOL    *
*        variable in GCTRAK common                                     *
*  IS    , set number                                                  *
*  ID    , detector number                                             *
*  IPA   , particle number                                             *
*  IT    , track number                                                *
*  DELOSS, energy loss in the step                                     *
*  PATH  , current track lenght                                        *
*  TRAPAR, track parameter, is the VECT vector in GCTRAK common        *
*  NUMVOL, integr array of numbers identifying the DETECTOR            *
*  IMEC  , integer indicating the LOSS mecanism position inside NMEC   *
*          vector, zero if not present                                 *
*                                                                      *
*                                                                      *
* Called by: GUSTEP                                                    *
* Author: Francesco Cafagna, 21/02/96 16.13.51                         *
*                                                                      *
************************************************************************
#include "gphit.inc"
#include "gcunit.inc"
#include "gptotr.inc"
#if defined(GPAMELA_GARFIELD)
#include "gctrak.inc"
#include "gckine.inc"
#endif
*
      INTEGER IACT,IS,ID,IPA,IT,NUMVOL(20),IMEC,IDHIT
      REAL DELOSS,PATH,TRAPAR(7),VHIT(NHTRD)
      LOGICAL SAVE
#if defined(GPAMELA_GARFIELD)
*
* Some variables need to call the GARFIELD interface
*
      CHARACTER*20 CHJUNK
      REAL GARXYZ(6),GARPAR(3),ECLTOT,RJUNK,RBJUNK(10)
      INTEGER IFAIL,NCLUS,IJUNK,IWBUFF
#endif
*
* Track is inside a volume
*
      IF(IACT.EQ.0) THEN
         IF(IMEC.GT.0) THEN
            SAVE = .TRUE.
            VHIT(7) = DELOSS + VHIT(7)
         ENDIF
      ENDIF
*
* Entering a new volume or is a new track
*
      IF(IACT.EQ.1) THEN
         SAVE = .FALSE.
         CALL VZERO(VHIT,NHTRD)
         DO I=1,3
            VHIT(I) = TRAPAR(I)
         ENDDO
         IF(IMEC.GT.0) THEN
            SAVE = .TRUE.
            VHIT(7) = DELOSS
         ENDIF
         VHIT(8) = PATH
         VHIT(9) = IPA
         VHIT(10)= TRAPAR(7)
         IDHIT = -1
      ENDIF
*
* Track is exiting current volume
*
      IF(IACT.EQ.2) THEN
         DO I=1,3
            VHIT(I+3) = TRAPAR(I)
#if defined(GPAMELA_INTER)
* interactive case, we need the average position
            VHIT(I) = (VHIT(I) + VHIT(I+3) )/2.
#endif
         ENDDO
#if defined(GPAMELA_GARFIELD)
*     Convert coordinates from the MARS to the local reference system
         IF(IMEC.GT.0) THEN
           CALL GMTOD(VHIT,GARXYZ,1)
           CALL GMTOD(TRAPAR,GARXYZ(4),1)
C     PRINT*,'GARXYZ,PATH=',(GARXYZ(I),I=1,6),PATH
*     CALL GPTOGAR(GARXYZ,GARXYZ)
*     CALL GPTOGAR(GARXYZ(4),GARXYZ(4))
*
*     Get the particle infos
*
           CALL GFPART(IPA,CHJUNK,IJUNK,GARPAR(1),GARPAR(2),
     +        RJUNK,RBJUNK,IWBUFF)
*
* Convert mass from GeV to Mev, calculate kinetic energy and convert it in MeV
*
           GARPAR(3) =
     +        (SQRT(TRAPAR(7)**2+GARPAR(1)**2)-GARPAR(1))*1000.
           GARPAR(1) = GARPAR(1)*1000.
*
*     Let's calculate the clusters
*
           CALL GPGACL(IPA,GARPAR,GARXYZ,ECLTOT,NCLUS)
*
*     fill the hit structure (energy in GeV)
*
           VHIT(13)=ECLTOT/1000.
C           VHIT(14)=FLOAT(NCLUS)
           VHIT(14)=INT(NCLUS)
C* ML: update the kinetic energy and then the momentum:
           GEKINT=GEKIN-VHIT(13)
           GEKIN=GEKINT
           GETOT=GEKIN +AMASS
           VECT(7)= SQRT((GETOT+AMASS)*GEKIN)
C  Routine to find bin number in kinetic energy table stored in ELOW(NEKBIN)
           CALL GEKBIN
C* END ML.
         ENDIF
#endif
         IF(IMEC.GT.0) THEN
            SAVE = .TRUE.
            VHIT(7) = DELOSS + VHIT(7)
         ENDIF
         VHIT(8) = PATH - VHIT(8)
*
*     Store the hit
*
         IF(IDHIT.EQ.-1) THEN
            IF(TROK) THEN
               VHIT(11)=ENPHTR*1.E-6
               VHIT(12)=NPHTR
            ELSE
               WRITE(CHMAIL,10000)
               CALL GMAIL(1,0)
            ENDIF
            IF(SAVE) CALL GSAHIT (IS,ID,IT,NUMVOL,VHIT,IDHIT)
            IF (IDHIT.EQ.0) THEN
               WRITE(CHMAIL,10100) IS,ID,IT
               CALL GMAIL(1,0)
            ENDIF
         ELSE
            WRITE(CHMAIL,10200) IS,ID,IT
            CALL GMAIL(1,0)
         ENDIF
      ENDIF
*
* Track is exiting the set up
*
      IF(IACT.EQ.3) THEN
*
* Check if the hit have been previously stored or not
*
         IF(IDHIT.LT.0) THEN
            DO I=1,3
               VHIT(I+3) = TRAPAR(I)
#if defined(GPAMELA_INTER)
* interactive case, we need the average position
               VHIT(I) = (VHIT(I) + VHIT(I+3) )/2.
#endif
            ENDDO
            IF(IMEC.GT.0) THEN
               SAVE = .TRUE.
               VHIT(7) = DELOSS + VHIT(7)
            ENDIF
            VHIT(8) = PATH - VHIT(8)
#if defined(GPAMELA_GARFIELD)
            IF(IMEC.GT.0) THEN
* Convert coordinates from the MARS to the local reference system
              CALL GMTOD(VHIT,GARXYZ,1)
              CALL GMTOD(TRAPAR,GARXYZ(4),1)
*     PRINT*,'GARXYZ,PATH=',(GARXYZ(I),I=1,6),PATH
*     CALL GPTOGAR(GARXYZ,GARXYZ)
*     CALL GPTOGAR(GARXYZ(4),GARXYZ(4))
*
*     Get the particle infos
*
              CALL GFPART(IPA,CHJUNK,IJUNK,GARPAR(1),GARPAR(2),
     +         RJUNK,RBJUNK,IWBUFF)
*
*  Convert mass from GeV to Mev, calculate kinetic energy and convert it in MeV
*
              GARPAR(3) =
     +         (SQRT(TRAPAR(7)**2+GARPAR(1)**2)-GARPAR(1))*1000.
              GARPAR(1) = GARPAR(1)*1000.
*
*     Let's calculate the clusters
*
              CALL GPGACL(IPA,GARPAR,GARXYZ,ECLTOT,NCLUS)
*
*     fill the hit structure (energy in GeV)
*
              VHIT(13)=ECLTOT/1000.
C              VHIT(14)=FLOAT(NCLUS)
              VHIT(14)=INT(NCLUS)
C ML: update the kinetic energy and then th momentum:
              GEKINT=GEKIN-VHIT(11)
              GEKIN=GEKINT
              GETOT=GEKIN +AMASS
              VECT(7)= SQRT((GETOT+AMASS)*GEKIN)
C Routine to find bin number in kinetic energy table stored in ELOW(NEKBIN)
              CALL GEKBIN
C END ML.
            ENDIF
#endif
*
* Store the hit
*
          IF(IDHIT.EQ.-1) THEN
            IF(TROK) THEN
               VHIT(11)=ENPHTR*1.E-6
               VHIT(12)=NPHTR
            ELSE
               WRITE(CHMAIL,10300)
               CALL GMAIL(1,0)
            ENDIF
          ENDIF
          IF(SAVE) CALL GSAHIT (IS,ID,IT,NUMVOL,VHIT,IDHIT)
          IF (IDHIT.EQ.0) THEN
            WRITE(CHMAIL,10400) IS,ID,IT
            CALL GMAIL(1,0)
          ENDIF
        ELSE
            WRITE(CHMAIL,10500) IS,ID,IT
            CALL GMAIL(1,0)
        ENDIF
      ENDIF

10000 FORMAT(' GPUTRD error: TR process and his tuning will be ignored')
10100 FORMAT(' GPUTRD error: HIT have not be stored for SET='
     +       ,I4,', DET=',I4,', TRACK=',I6)
10200 FORMAT(' GPUTRD error: HIT have not be initialized for SET='
     +       ,I4,', DET=',I4,', TRACK=',I6)
10300 FORMAT(' GPUTRD error: TR process and his tuning will be ignored')
10400 FORMAT(' GPUTRD error: HIT have not be stored for SET='
     +       ,I4,', DET=',I4,', TRACK=',I6)
10500 FORMAT(' GPUTRD error: HIT have not be initialized for SET='
     +       ,I4,', DET=',I4,', TRACK=',I6)
      RETURN
      END
