*
* $Id$
*
* $Log$
*
*CMZ :  3.00/00 04/04/2001  15.19.56  by  Marialuigia Ambriola
*CMZ :  2.03/00 06/11/2000  01.58.14  by  Francesco Cafagna
*CMZ :  2.02/00 12/10/2000  18.37.54  by  Francesco Cafagna
*CMZU:  2.01/00 04/04/2000  17.20.52  by  Unknown
*CMZ :  2.00/00 01/03/2000  21.21.39  by  Francesco Cafagna
*CMZ :  1.02/00 31/12/99  17.53.56  by  Francesco Cafagna
*CMZ :  1.01/00 04/05/96  09.50.23  by  Francesco Cafagna
*CMZ :  1.00/03 17/04/96  17.35.58  by  Francesco Cafagna
*CMZ :  1.00/02 05/04/96  09.27.01  by  Francesco Cafagna
*-- Author :    Francesco Cafagna   05/12/95
      SUBROUTINE GPDAT
************************************************************************
*                                                                      *
* Data initialization for PAMELA.                                      *
* Geometry parameter and general constant are definited                *
* Called by: UGINIT                                                    *
* Author: Francesco Cafagna, 05/12/95 13.21.27                         *
*                                                                      *
************************************************************************
#include "gclist.inc"
#include "gcflag.inc"
#include "gcunit.inc"
*
#include "gpgeo.inc"
#include "gpmed.inc"
#include "gpmat.inc"
#include "gpkey.inc"
#include "gphit.inc"
#include "gpdigi.inc"
#include "gpfield.inc"
#include "gprzfi.inc"
#include "gpunit.inc"
#include "gpal.inc"
#include "gpg10.inc"
#include "gpn2g.inc"
#include "gpw.inc"
#include "gpfe.inc"
#include "gpkap.inc"
#include "gptrad.inc"
#include "gproa.inc"
#include "gpcp.inc"
#include "gpscint.inc"
#include "gpsica.inc"
#include "gpsitr.inc"
#include "gpxe.inc"
#include "gpaer.inc"
#include "gpgene.inc"
*EM:
#include "gpce.inc"
#include "gpw2.inc"
#include "gpg10c.inc"
*END: EM.
*ML:
#include "gpaltr.inc"
*END ML.
*
      INTEGER ICP
      LOGICAL TEST,GPLOOK
      CHARACTER*4 CJUNK
*
* Data statements
*
#include "gpdgeo.inc"
#include "gpdmat.inc"
#include "gpdhit.inc"
#include "gpddigi.inc"
#include "gpdal.inc"
#include "gpdg10.inc"
#include "gpdn2g.inc"
#include "gpdw.inc"
#include "gpdfe.inc"
#include "gpdkap.inc"
#include "gpdtrad.inc"
#include "gpdroa.inc"
#include "gpdcp.inc"
#include "gpdscin.inc"
#include "gpdsica.inc"
#include "gpdsitr.inc"
#include "gpdxe.inc"
#include "gpdaer.inc"
*EM:
#include "gpdce.inc"
#include "gpdw2.inc"
#include "gpdg10c.inc"
*END: EM.
*ML:
C
C Read the attenuation lengths for the TRD from the file
C
      READ(LUTR) ENATT,ATTTRD
      CALL VZERO(EY,115)
*END ML.
C
C Set quantities from DATA cards
C
C
C Number of plane without tungsten
C
      ICP = -1
      NCANOW = 1
      DO I=1,3
         CALL UHTOC(NCAL(I),4,CJUNK,4)
         ICP = ICFIND('NOW',CJUNK,1,4)
         IF(ICP.EQ.1) NCANOW = ICDECI(CJUNK,4,4)
      ENDDO
C
C CWN booking and filling
C
      IF(HBOOK.EQ.1) THEN
         HBK = .TRUE.
      ELSE
         HBK = .FALSE.
      ENDIF
C
c Look for longer CWN-tplue selection
C
      IF(HLON.EQ.1) THEN
         HLONG = .TRUE.
      ELSE
         HLONG =.FALSE.
      ENDIF

C
C Volume to be subtracted via NDET card
C
      SPH=.NOT.GPLOOK('SPHE',NDET,8)
C
      TOF=.NOT.GPLOOK('TOF ',NDET,8)
      HTOF = .FALSE.
      IF(HBK.AND.TOF) HTOF = .TRUE.
C
      TRD=.NOT.GPLOOK('TRD ',NDET,8)
      HTRD = .FALSE.
      IF(HBK.AND.TRD) HTRD = .TRUE.
C
      CAT =.NOT.GPLOOK('CAT ',NDET,8)
      HCAT = .FALSE.
      IF(HBK.AND.CAT) HCAT = .TRUE.
C
      CAS=.NOT.GPLOOK('CAS ',NDET,8)
      HCAS = .FALSE.
      IF(HBK.AND.CAS) HCAS = .TRUE.
C
      SPE=.NOT.GPLOOK('SPE ',NDET,8)
      HSPE = .FALSE.
      IF(HBK.AND.SPE) HSPE = .TRUE.
C
      CAL=.NOT.GPLOOK('CAL ',NDET,8)
      HCAL = .FALSE.
      IF(HBK.AND.CAL) HCAL = .TRUE.
C
C E-ntple or physics disablign for single detectors
C
      IF (CAS) THEN
         CAS=.NOT.GPLOOK('NVOL',NDECAS,3)
         PCAS=.NOT.GPLOOK('NPHY',NDECAS,3)
         HCAS=.NOT.GPLOOK('NHBK',NDECAS,3)
         IF(.NOT.CAS) PCAS = .FALSE.
         IF(.NOT.PCAS)HCAS = .FALSE.
         IF(.NOT.HBK) HCAS = .FALSE.
      ELSE
         PCAS = .FALSE.
         HCAS = .FALSE.
      ENDIF
      IF (CAT) THEN
         CAT=.NOT.GPLOOK('NVOL',NDECAT,3)
         PCAT=.NOT.GPLOOK('NPHY',NDECAT,3)
         HCAT=.NOT.GPLOOK('NHBK',NDECAT,3)
         IF(.NOT.CAT) PCAT = .FALSE.
         IF(.NOT.PCAT)HCAT = .FALSE.
         IF(.NOT.HBK) HCAT = .FALSE.
      ELSE
         PCAT = .FALSE.
         HCAT = .FALSE.
      ENDIF
      IF (TOF) THEN
         TOF=.NOT.GPLOOK('NVOL',NDETOF,3)
         PTOF=.NOT.GPLOOK('NPHY',NDETOF,3)
         HTOF=.NOT.GPLOOK('NHBK',NDETOF,3)
         IF(.NOT.TOF) PTOF = .FALSE.
         IF(.NOT.PTOF)HTOF = .FALSE.
         IF(.NOT.HBK) HTOF = .FALSE.
      ELSE
         PTOF = .FALSE.
         HTOF = .FALSE.
      ENDIF
      IF (TRD) THEN
         TRD=.NOT.GPLOOK('NVOL',NDETRD,3)
         PTRD=.NOT.GPLOOK('NPHY',NDETRD,3)
         HTRD=.NOT.GPLOOK('NHBK',NDETRD,3)
         IF(.NOT.TRD) PTRD = .FALSE.
         IF(.NOT.PTRD)HTRD = .FALSE.
         IF(.NOT.HBK) HTRD = .FALSE.
      ELSE
         PTRD = .FALSE.
         HTRD = .FALSE.
      ENDIF
      IF (SPE) THEN
         SPE=.NOT.GPLOOK('NVOL',NDESPE,3)
         PSPE=.NOT.GPLOOK('NPHY',NDESPE,3)
         HSPE=.NOT.GPLOOK('NHBK',NDESPE,3)
         IF(.NOT.SPE) PSPE = .FALSE.
         IF(.NOT.PSPE)HSPE = .FALSE.
         IF(.NOT.HBK) HSPE = .FALSE.
      ELSE
         PSPE = .FALSE.
         HSPE = .FALSE.
      ENDIF
      IF (CAL) THEN
         CAL=.NOT.GPLOOK('NVOL',NDECAL,3)
         PCAL=.NOT.GPLOOK('NPHY',NDECAL,3)
         HCAL=.NOT.GPLOOK('NHBK',NDECAL,3)
         IF(.NOT.CAL) PCAL = .FALSE.
         IF(.NOT.PCAL)HCAL = .FALSE.
         IF(.NOT.HBK) HCAL = .FALSE.
      ELSE
         PCAL = .FALSE.
         HCAL = .FALSE.
      ENDIF
      IF (TS4) THEN
         TS4=.NOT.GPLOOK('NVOL',NDES4,3)
         PS4=.NOT.GPLOOK('NPHY',NDES4,3)
         HS4=.NOT.GPLOOK('NHBK',NDES4,3)
         IF(.NOT.TS4)  PS4 = .FALSE.
         IF(.NOT.PS4) HS4 = .FALSE.
         IF(.NOT.HBK) HS4 = .FALSE.
      ELSE
         PS4 = .FALSE.
         HS4 = .FALSE.
      ENDIF
*
* Look for special tracking media settings
*
      TEST =GPLOOK('AL  ',NSPTM,17)
      IF(TEST) CALL GPAL
      TEST =GPLOOK('G10 ',NSPTM,17)
      IF(TEST)  CALL GPG10
*EM:
      TEST =GPLOOK('G10C',NSPTM,17)
      IF(TEST)  CALL GPG10C
*END: EM.
      TEST =GPLOOK('N2G ',NSPTM,17)
      IF(TEST)  CALL GPN2G
      TEST =GPLOOK('W   ',NSPTM,17)
      IF(TEST)  CALL GPW
*EM:
      TEST =GPLOOK('W2  ',NSPTM,17)
      IF(TEST)  CALL GPW2
*END: EM.
      TEST =GPLOOK('FE  ',NSPTM,17)
      IF(TEST)  CALL GPFE
      TEST =GPLOOK('KAP ',NSPTM,17)
      IF(TEST)  CALL GPKAP
      TEST =GPLOOK('TRAD',NSPTM,17)
      IF(TEST)  CALL GPTRAD
      TEST =GPLOOK('ROA ',NSPTM,17)
      IF(TEST)  CALL GPROA
      TEST =GPLOOK('CP  ',NSPTM,17)
      IF(TEST)  CALL GPCP
      TEST =GPLOOK('SCIN',NSPTM,17)
      IF(TEST)  CALL GPSCINT
      TEST =GPLOOK('SICA',NSPTM,17)
      IF(TEST)  CALL GPSICA
      TEST =GPLOOK('SITR',NSPTM,17)
      IF(TEST)  CALL GPSITR
      TEST =GPLOOK('XE  ',NSPTM,17)
      IF(TEST)  CALL GPXE
*EM:
      TEST =GPLOOK('CE  ',NSPTM,17)
      IF(TEST)  CALL GPCE
*END EM.
*
* Look for parameter to be always atomaticaly calculated by GEANT
*
      ATMAXFD = GPLOOK('TMAX',NAUTMED,5)
      ASTEMAX = GPLOOK('STMA',NAUTMED,5)
      ADEEMAX = GPLOOK('DEEM',NAUTMED,5)
      AEPSIL  = GPLOOK('EPSI',NAUTMED,5)
      ASTMIN  = GPLOOK('STMI',NAUTMED,5)
*
* Look for hadronic package to be used
*
      FLUKA =GPLOOK('FLUK',NHPAK,1)
*
* Look for Gaussian straggling parameter
* Please note that is updated only if ISSTR is one !
*
      IF(ISSTR.EQ.1) THEN
         IF(SSTR(1).NE.0) D2SICA = SSTR(1)
         IF(SSTR(2).NE.0) D2SITR = SSTR(2)
      ENDIF
*
* Look for calorimeter planes to be excluded
*
      DO I=1,2*NCAPL
         IF(NCPL(I).GT.0) THEN
            IF(NCPL(I).LE.2*NCAPL) THEN
               NOCAPL(NCPL(I))=.TRUE.
            ELSE
               WRITE(CHMAIL,10100) NCPL(I)
               CALL GMAIL(1,0)
            ENDIF
         ENDIF
      ENDDO
*
* Look for calorimeter silicon detectors to be excluded
*
      DO I=1,9
         IF(NCSI(I).GT.0) THEN
            IF(NCSI(I).LE.9) THEN
               NOCASI(NCSI(I))=.TRUE.
            ELSE
               WRITE(CHMAIL,10000) NCSI(I)
               CALL GMAIL(1,0)
            ENDIF
         ENDIF
      ENDDO
*
* Look for calorimeter absorber planes to be excluded
*
      DO I=1,NCAPL
         IF(NCAB(I).GT.0) THEN
            IF(NCAB(I).LE.22) THEN
               NOCAAB(NCAB(I))=.TRUE.
            ELSE
               WRITE(CHMAIL,10200) NCAB(I)
               CALL GMAIL(1,0)
            ENDIF
         ENDIF
      ENDDO
*
* Origin of PAMELA in the MARS GEANT coordinates system
*
      XYZPOR(1) = 0.
      XYZPOR(2) = 0.
      XYZPOR(3) = 0.
C
C TOF dimensions
C
      S11X(1)=S11(1)/NS11X
      S11X(2)=S11(2)
      S11X(3)=S11(3)-2.*SCTIC
      SC11(1)=S11(1)
      SC11(2)=S11(2)
      SC11(3)=SCTIC
      S12Y(1)=S12(1)
      S12Y(2)=S12(2)/NS12Y
      S12Y(3)=S12(3)-2.*SCTIC
      SC12(1)=S12(1)
      SC12(2)=S12(2)
      SC12(3)=SCTIC
      S21X(1)=S21(1)/NS21X
      S21X(2)=S21(2)
      S21X(3)=S21(3)-2.*SCTIC
      SC21(1)=S21(1)
      SC21(2)=S21(2)
      SC21(3)=SCTIC
      S22Y(1)=S22(1)
      S22Y(2)=S22(2)/NS22Y
      S22Y(3)=S22(3)-2.*SCTIC
      SC22(1)=S22(1)
      SC22(2)=S22(2)
      SC22(3)=SCTIC
      S31Y(1)=S31(1)
      S31Y(2)=S31(2)/NS31Y
      S31Y(3)=S31(3)-2.*SCTIC
      SC31(1)=S31(1)
      SC31(2)=S31(2)
      SC31(3)=SCTIC
      S32X(1)=S32(1)/NS32X
      S32X(2)=S32(2)
      S32X(3)=S32(3)-2.*SCTIC
      SC32(1)=S32(1)
      SC32(2)=S32(2)
      SC32(3)=SCTIC
      ZS4 = ZCAL - CALS(3) - S4(3)
*
* Calculate TRD straw box dimentions add some dead space
*
C #       TRBS(1) = (16.+0.5)*TRSO(2)
C #       TRBS(2) = TRSO(3)
C # +SELF,IF=UNIX
C #       TRBS(3) = 0.246+(1+COS(30./180.*ACOS(-1.)))*TRSO(2)
C # +SELF.
C #       TRBS(3) = 0.246+(1+COSD(30.))*TRSO(2)
C # +SELF.
*
* Calculate TRD radiator thickness
*
      TRRA(3) = (TRBS(3) - (TRSO(2)*( 1 + COS(30./180.*ACOS(-1.))))) +
     +           TRFR(3)
*
* Calculate TRD virtual box dimensions
*
      TRDB(1) = TRFR(1)
      TRDB(2) = TRFR(2)
*      TRDB(3) = 9*TRBS(3) + 9*TRFR(3) + 2* TRAN(3) + TRDT(3)
      TRDB(3) = 10*TRFR(3) + 11* TRAN(3) + TRDT(3)
*
* Calculate the silicon ladder and spectrometer plane dimensions
*
      TRSL(1) = TSPA(1) + 2.*TRCP(1)
      TRSL(3) = TRCP(3)
      TRPL(1) = TRSL(1)*NTRSL
*
* Calculate spectrometer box Z dimentions
*
      SPEB(3) = 5.*MGPL(3) + 6.*TRPL(3)
*
* Calculate the calorimeter planes dimentions, both for actual and dummy planes
*
C fc mod by caf on 20-10-1999
C      CAPL(3) = ((2.*CALB(3)/NCAPL - 2*CAAB(3))/2.)/2.
C end of mod
*EM:
c      CAPL(3) = (2.*CG10(3) + 2*CASI(3) + DCASIZ/2. )/2.
      CAPL(3) = (DCASIZ/2. + 2*CASI(3) + 2.*CAKP(3) +
     + 2.*CAGL(3)+ 2.*C10C(3) + 2*CAKA(3))/2.
      CAPD(3) = CAPL(3)
C
C We have repectively: N2, Si, Kapton, Glue, G10C, Kaolinite
C
*END EM.
      CASIOFF = (2.*CASI(1) - CALSTWID*NCASTR)/2. - CASI(1)
*
* Calculate the inner shell dimentions and the TOP shell maximum theta angle
* The inner shel is two centimeters from the external and 1 mm thick
*
      SHEI(2) = SHEL(2) - 2.
      SHEI(1) = SHEI(2) - .1
      TSPH(1) = TSPH(2) - .25
#if defined(GPAMELA_UNIX)
      TSPH(4) = ASIN(SHEL(2)/TSPH(2)) * 180./ACOS(-1.)
#endif
#if !defined(GPAMELA_UNIX)
      TSPH(4) = ASIND(SHEL(2)/TSPH(2))
#endif
*
* Quote definition for the main volumes
*
* Note that the quote are in PAMELA System of reference. They are transformed
* in GEANT coordinates system
*
*EM:
C      ZCAS = ZSPEC
*END EM.
*
* The shells heights are the same of the TRD height
*
      SHEL(3) = (ZTRD + TRDB(3) +PAME(3) )/2.
*
* Add user defined hieght to the external shell
*
      SHEL(3) = SHEL(3) + DZSH
      SHEI(3) = SHEL(3)
      ZSHEL = SHEL(3) - PAME(3)
      ZSHEI = SHEI(3) - PAME(3)
      ZTFLA = ZSHEL + SHEL(3) + TFLA(3)
#if defined(GPAMELA_UNIX)
      ZTSPH = (2.*SHEI(3)-PAME(3))-TSPH(2)*COS(TSPH(4)/180.*ACOS(-1.))
#endif
#if !defined(GPAMELA_UNIX)
      ZTSPH = (2.*SHEI(3)-PAME(3)) - TSPH(2) * COSD(TSPH(4))
#endif
*
* Origin of the spectrometer in the MARS GEANT coordinates system
*
      XYZSOR(1) = 0.
      XYZSOR(2) = 0.
      XYZSOR(3) = ZSPEC
*
* Calculate generation surface & quote based on the S11 dimentions
*
      IF(XYZGEN(1).NE.-1111) THEN
         XGEN = XYZGEN(1)
      ELSE
         XGEN = -S11(1)
      ENDIF
      IF(XYZGEN(2).NE.-1111) THEN
         YGEN = XYZGEN(2)
      ELSE
         YGEN = -S11(2)
      ENDIF
      IF(XYZGEN(3).NE.-1111) THEN
         ZGEN = XYZGEN(3)
      ELSE
         ZGEN = ZS11 + 2.
      ENDIF
      IF(XYZGEN(4).NE.-1111) THEN
         XDGEN = XYZGEN(4)
      ELSE
         XDGEN = 2.*S11(1)
      ENDIF
      IF(XYZGEN(5).NE.-1111) THEN
         YDGEN = XYZGEN(5)
      ELSE
         YDGEN= 2.*S11(2)
      ENDIF
C
C Set max and min generation angles
C
      THMIN = THGEN(1)
      THMAX = THGEN(2)
      PHMIN = PHGEN(1)
      PHMAX = PHGEN(2)
C
C Tracking media parameters. These have been set following M. Boezio
C suggestions, to be optimized for silicon detectors simulation
*ML:
C The EPSILGP variable will be set in GPMED routine.
* The EPSILGP variable is set now and, if it is the case, it is modified
* in GPMED routine.
*ML.
C The FFIELD variable is passed by a DATA CARD
C
      FIELDMGP = 40.
      IF(VTMAXFD.GT.0.) THEN
         TMAXFDGP = VTMAXFD
*ML:
         BTMAXFD = .TRUE.
*END ML.
      ELSE
         TMAXFDGP = 1.
      ENDIF
      IF(VSTEMAX.GT.0.) THEN
         STEMAXGP = VSTEMAX
*ML:
         BSTEMAX = .TRUE.
*END ML.
      ELSE
         STEMAXGP = 0.05
      ENDIF
      IF(VDEEMAX.GT.0.) THEN
         DEEMAXGP = VDEEMAX
*ML:
         BDEEMAX = .TRUE.
*END ML.
      ELSE
         DEEMAXGP = 0.05
      ENDIF
      IF(VSTMIN.GT.0.) THEN
         STMINGP = VSTMIN
*ML:
         BSTMIN = .TRUE.
*END ML.
      ELSE
         STMINGP = 0.0005
      ENDIF
*ML:
      IF(VEPSIL.GT.0.) THEN
         EPSILGP = VEPSIL
         BEPSIL = .TRUE.
      ENDIF
*END ML.
C
C Read magnetic field map
C
      IF (FFIELD.EQ.1) THEN
C fc          DO IX=1,41
C fc             DO IY=1,41
C fc                DO IZ=1,31
C fc                   READ(LUFMAP,ERR=10) FX(IX,IY,IZ),FY(IX,IY,IZ),
C fc      +                                FZ(IX,IY,IZ)
C fc                END DO
C fc             END DO
C fc          END DO
C fc    10    PRINT*,'ERROR READING MAP FILE',IX,IY,IZ
         CALL HCDIR('//FIELD',' ')
         CALL HRIN(30,9999,0)
         CALL HPRNT(30)
         CALL HBNAME(30,' ',0,'$CLEAR')
         CALL HBNAME(30,'FIELD',F,'$SET')
         CALL HNOENT(30,IPTOT)
         IP=0
         DO IX=1,41
            DO IY=1,41
               DO IZ=1,31
                  IP=IP+1
                  CALL HGNT(30,IP,IERR)
                  FX(IX,IY,IZ)=F(1)
                  FY(IX,IY,IZ)=F(2)
                  FZ(IX,IY,IZ)=F(3)
               END DO
            END DO
         END DO
         CALL HREND('FIELD')
      ENDIF
*
* Set run and first event number from RUNG GEANT data card
*
      IRUN = IDRUN
      IEVNT = IDEVT
C
C User action
C
      CALL GPUDAT
10000 FORMAT(' GPDDAT: NCPL out of range, NCPL=',I8,'. Not used')
10100 FORMAT(' GPDDAT: NCSI out of range, NCSI=',I8,'. Not used')
10200 FORMAT(' GPDDAT: NCAB out of range, NCAB=',I8,'. Not used')
      RETURN
      END
