*
* $Id: gpdat.F,v 3.22 2006/11/30 12:28:39 cafagna Exp $
*
* $Log: gpdat.F,v $
* Revision 3.22  2006/11/30 12:28:39  cafagna
* CARD geometry updated. New S1 "table", legs included, added. Top magnetic screen has been added as well.
*
* Revision 3.21  2006/11/28 10:26:14  pam-ba
* S3 positioning completed
*
* Revision 3.20  2006/11/16 18:45:28  pam-ba
* Simulated an aluminum container for S4
*
* Revision 3.19  2006/11/16 12:01:12  pam-ba
* S4 dimensions corrected.
*
* Revision 3.18  2006/11/10 11:39:33  pam-ba
* S2 and S1 z-positions corrected, He3 and plystyrene mixture added, Top Plate geometry simulated and titanium mixture added.
*
* Revision 3.17  2006/10/13 16:36:58  pam-ba
* Added a new material, the cadmium, for ND
*
* Revision 3.16  2006/10/12 11:11:20  pam-ba
* ND geometry updated.
*
* Revision 3.15  2006/06/30 15:38:15  pam-ba
* S22 and S12 heights positioned in GPAMELA at the nominal heights in PAMELA (see document: Main geometrical parameters of the PAMELA sub-detectors, 20 December 2005)
*
* Revision 3.14  2006/06/06 10:14:38  pam-ba
* gpmyl.inc and gpdmyl.inc included in gpdat.F, polystyrene data card presence test added in gpdat.F and polystyrene number assigned in gpmat.F
*
* Revision 3.13  2006/06/05 13:56:16  pamela
* Gigantic resonance added for gamma enetering in the calorimeter absorber
*
* Revision 3.12  2006/05/18 10:52:31  pam-ba
* TOF geometry completed and a new material, the polystyrene (density 35 g/l), added
*
* Revision 3.11  2006/05/02 12:24:43  bottai
*  clean a global variable for strip simulation
*
* Revision 3.10  2006/05/02 11:46:13  pam-ba
* TOF geometry and position updated and a new material, the mylar, added
*
* Revision 3.9  2006/04/10 11:07:42  cafagna
* GEN data card updated, ZDGEN added
*
* Revision 3.8  2006/04/09 23:28:32  cafagna
* Several new things, among this: ND and CARD
*
* Revision 3.7  2005/12/20 12:21:03  cafagna
* gpnd directory added along with ND files
*
* Revision 3.6  2005/12/14 03:12:39  cafagna
* Neutron detector added. Geometry and GPCALOR package
*
* Revision 3.5  2005/12/05 12:15:20  pam-ba
* new spectrometer geometry and internal magnetic field
*
* Revision 3.4  2005/07/25 11:53:20  cafagna
* Several updates. See history for details
*
* Revision 3.3  2005/06/21 02:41:43  cafagna
* Major modification to the geometry and to the random number chain
*
* Revision 3.2  2002/12/04 15:12:52  pamela
* CAS and CALO geometry modified
*
* Revision 3.1.1.1  2002/07/11 16:02:00  cafagna
* First GPAMELA release on CVS
*
*
*               14/04/2005  13:38     by  Marialuigia Ambriola
*               16/07/2004  16.44.06  by  Alexey Bakaldin
*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 "gprspe.inc"
#include "gpdrspe.inc"
#include "gpgeo.inc"
#include "gpmed.inc"
#include "gpmat.inc"
#include "gpkey.inc"
#include "gphit.inc"
#include "gpdigi.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.
#include "gpmyl.inc"
#include "gppol.inc"
#include "gpcad.inc"
#include "gpch2.inc"
#include "gphe3.inc"
#include "gptav.inc"
*ML:
#include "gpaltr.inc"
#include "gpseed.inc"
*END ML.
#include "gpstripspe.inc"
#include "gpgig.inc"
#include "gppgmec.inc"
#include "gctrak.inc"

      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.
#include "gpdmyl.inc"
#include "gpdpol.inc"
#include "gpdcad.inc"
#include "gpdch2.inc"
#include "gpdhe3.inc"
#include "gpdtav.inc"
*ML:
C
C Read the attenuation lengths for the TRD from the file
C
      READ(LUTR) ENATT,ATTTRD
      CLOSE (LUTR)
      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
*ML:29/4/04
C
C Look for the input seeds of the random number generator 
C
C # Please note this superseeds the GEANT defaults set by RNDM card
C # or the ones set by GEANT in ginit.F
C # HFSF = 1, then just read the seeds used to initialize a previous run, 
C #           i.e. repeat the random chain
      IF(HSEED.EQ.1) THEN
         READ(LUFSF) INSEED
         NRNDM(1)=INSEED(1)
         NRNDM(2)=INSEED(2)
         CALL GRNDMQ(NRNDM(1),NRNDM(2),1,'S')
C # HFSF = 2, then just read the last seeds used in a previous run, 
C #           i.e. continue the random chain
      ELSEIF(HSEED.EQ.2) THEN
         READ(LUFSF) INSEED
         NRNDM(1)=INSEED(3)
         NRNDM(2)=INSEED(4)
         CALL GRNDMQ(NRNDM(1),NRNDM(2),1,'S')
      ELSE
C # HFSF not set, use the default GEANT setting
         CALL GRNDMQ(NRNDM(1),NRNDM(2),0,'G')
      ENDIF
      OUTSEED(1)=NRNDM(1)
      OUTSEED(2)=NRNDM(2)
*END ML.

C
C Volume to be subtracted via NDET card
C
      SPH=.NOT.GPLOOK('SPHE',NDET,10)
C
      TOF=.NOT.GPLOOK('TOF ',NDET,10)
      HTOF = .FALSE.
      IF(HBK.AND.TOF) HTOF = .TRUE.
C
      TRD=.NOT.GPLOOK('TRD ',NDET,10)
      HTRD = .FALSE.
      IF(HBK.AND.TRD) HTRD = .TRUE.
C
      CAT =.NOT.GPLOOK('CAT ',NDET,10)
      HCAT = .FALSE.
      IF(HBK.AND.CAT) HCAT = .TRUE.
C
      CAS=.NOT.GPLOOK('CAS ',NDET,10)
      HCAS = .FALSE.
      IF(HBK.AND.CAS) HCAS = .TRUE.
C
      CARD=.NOT.GPLOOK('CARD',NDET,10)
      HCARD = .FALSE.
      IF(HBK.AND.CARD) HCARD = .TRUE.
C
      SPE=.NOT.GPLOOK('SPE ',NDET,10)
      HSPE = .FALSE.
      IF(HBK.AND.SPE) HSPE = .TRUE.
C
      CAL=.NOT.GPLOOK('CAL ',NDET,10)
      HCAL = .FALSE.
      IF(HBK.AND.CAL) HCAL = .TRUE.
C
      ND=.NOT.GPLOOK('ND  ',NDET,10)
      HND = .FALSE.
      IF(HBK.AND.ND) HND = .TRUE.
c$$$
c$$$ TRD AND CARD cannot be activate at the same time. 
c$$$ Precedence to the CARD is given. 
c$$$
      IF(CARD.AND.TRD) THEN 
         TRD = .FALSE.
         HTRD = .FALSE.
      ENDIF
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 (CARD) THEN
         CARD=.NOT.GPLOOK('NVOL',NDECARD,3)
         PCARD=.NOT.GPLOOK('NPHY',NDECARD,3)
         HCARD=.NOT.GPLOOK('NHBK',NDECARD,3)
         IF(.NOT.CARD) PCARD = .FALSE.
         IF(.NOT.PCARD)HCARD = .FALSE.
         IF(.NOT.HBK)  HCARD = .FALSE.
      ELSE
         PCARD = .FALSE.
         HCARD = .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
      IF (ND) THEN
         ND=.NOT.GPLOOK('NVOL',NDEND,3)
         PND=.NOT.GPLOOK('NPHY',NDEND,3)
         HND=.NOT.GPLOOK('NHBK',NDEND,3)
         IF(.NOT.ND)  PND = .FALSE.
         IF(.NOT.PND) HND = .FALSE.
         IF(.NOT.HBK) HND = .FALSE.
      ELSE
         PND = .FALSE.
         HND = .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
      TEST =GPLOOK('POL ',NSPTM,17)
      IF(TEST) CALL GPPOL
      TEST =GPLOOK('MYL ',NSPTM,17)
      IF(TEST) CALL GPMYL
      TEST =GPLOOK('CAD ',NSPTM,17)
      IF(TEST) CALL GPCAD
*EM:
      TEST =GPLOOK('CE  ',NSPTM,17)
      IF(TEST)  CALL GPCE
*END EM.
      TEST =GPLOOK('CH2 ',NSPTM,17)
      IF(TEST) CALL GPCH2
      TEST =GPLOOK('HE3 ',NSPTM,17)
      IF(TEST) CALL GPHE3
      TEST =GPLOOK('TAV ',NSPTM,17)
      IF(TEST) CALL GPTAV
*
* 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)
      GPCALOR =GPLOOK('GCAL',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
c      S11M(1)=S11Y(1)*NS11Y+2.*(HMYLTOF/2.)
      S11M(1)=S11Y(1)+2.*(HMYLTOF/2.)
      S11M(2)=S11Y(2)
      S11M(3)=S11Y(3)+2.*(HMYLTOF/2.)
      S12M(1)=S12X(1)
c      S12M(2)=S12X(2)*NS12X+2.*(HMYLTOF/2.)
      S12M(2)=S12X(2)+2.*(HMYLTOF/2.)
      S12M(3)=S12X(3)+2.*(HMYLTOF/2.)
      S21M(1)=S21X(1)
c      S21M(2)=S21X(2)*NS21X+2.*(HMYLTOF/2.)
      S21M(2)=S21X(2)+2.*(HMYLTOF/2.)
      S21M(3)=S21X(3)+2.*(HMYLTOF/2.)
c      S22M(1)=S22Y(1)*NS22Y+2.*(HMYLTOF/2.)
      S22M(1)=S22Y(1)+2.*(HMYLTOF/2.)
      S22M(2)=S22Y(2)
      S22M(3)=S22Y(3)+2.*(HMYLTOF/2.)
C      S31M(1)=S31Y(1)*NS31Y+3.*(HMYLTOF/2.)
      S31M(1)=S31Y(1)+3.*(HMYLTOF/2.)
      S31M(2)=S31Y(2)
      S31M(3)=S31Y(3)+3.*(HMYLTOF/2.)
      S32M(1)=S32X(1)
C      S32M(2)=S32X(2)*NS32X+3.*(HMYLTOF/2.)
      S32M(2)=S32X(2)+3.*(HMYLTOF/2.)
      S32M(3)=S32X(3)+3.*(HMYLTOF/2.)
*
c      S1A(1)=S12M(1)
c      S1A(2)=S12M(2)
      S1A(1)=S11M(1)*NS11Y
      S1A(2)=S12M(2)*NS12X
      S1A(3)=2.*(HGLUTOF/2.)+HPVCMIN/2.+S12M(3)+S11M(3)+HPVCMAX/2.
      POL1(1)=S1A(1)
      POL1(2)=S1A(2)
      POL1(3)=HPVCMAX/2.
      POL2(1)=S1A(1)
      POL2(2)=S1A(2)
      POL2(3)=HPVCMIN/2.
      S1(1)=S1A(1)
      S1(2)=S1A(2)
      S1(3)=2.*(HALTOF/2)+S1A(3)
*
c      S2A(1)=S21M(1)
c      S2A(2)=S22M(2)
      S2A(1)=S22M(1)*NS22Y
      S2A(2)=S21M(2)*NS21X
      S2A(3)=2*(HGLUTOF/2)+S22M(3)+S21M(3)+HPVCMIN/2
      POLY(1)=S2A(1)
      POLY(2)=S2A(2)
      POLY(3)=HPVCMIN/2.
      S2(1)=S2A(1)
      S2(2)=S2A(2)
      S2(3)=2*(HALTOF/2)+S2A(3)
*
C      S3(1)=S32M(1)
C      S3(2)=S31M(2)
      S3(1)=S31M(1)*NS31Y
      S3(2)=S32M(2)*NS32X
      S3(3)=HAIRS3/2.+S31M(3)+S32M(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
*
c      TRRA(3) = (TRBS(3) - (TRSO(2)*( 1 + COS(30./180.*ACOS(-1.))))) +
c     +           TRFR(3)
*
* Calculate TRD virtual box dimensions
*
c ml:17/11/04:
c      TRDB(1) = TRFR(1)
c      TRDB(2) = TRFR(2)
      TRDB(1) = TRR0(1)
      TRDB(2) = TRR0(2)
c end ml.
*      TRDB(3) = 9*TRBS(3) + 9*TRFR(3) + 2* TRAN(3) + TRDT(3)
c ml: 11/11/04:
c      TRDB(3) = 10*TRFR(3) + 11* TRAN(3) + TRDT(3)
c      TRDB(3)=TRR0(3)+9*TRAN(3)+8*TRFR(3)+TRR2(3)
      TRDB(3)=(TRR0(3)-0.05)+9*TRAN(3)+8*TRFR(3)+TRFM(3)+TRFL(3)
c end ml.
*
* Calculate the silicon ladder and spectrometer plane dimensions
*
c ml: 26/11/04:
c      TRSL(1) = TSPA(1) + 2.*TRCP(1)
c      TRSL(3) = TRCP(3)
c      TRPL(1) = TRSL(1)*NTRSL
c ml: 27/12/04:
c      TRPL(1) = TRSL(1)*NTRSL+4*TRCP(1)
c end ml.
c end ml.
*
* Calculate spectrometer box Z dimentions
*
cc 20/04/05:
c      SPEB(3)=5.*MGFR(3)+TRPB(3)
CC end ml.
*
* 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.
*F.V.&ML:
C      CAPL(3) = (DCASIZ/2. + 2*CASI(3) + 2.*CAKP(3) +
C     + 2.*CAGL(3)+ 2.*C10C(3) + 2*CAKA(3))/2.
C      CAPD(3) = CAPL(3)
* MA
C #       CAPL(3) =  (DCASIZ/2. + 2*CASI(3))/2.
      CAPL(3) =  CASI(3)
* end MA
      CANS(3) = CAKP(3) + CAGL(3)+ C10C(3) + CAKA(3)
      CAPD(3) = CANS(3)
*END F.V.&ML.
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
*

*AB:
C      SHEI(2) = SHEL(2) - 2.
      SHEI(2)=0.
C      SHEI(1) = SHEI(2) - .1
      SHEI(1)=0.
C      TSPH(1) = TSPH(2) - .25
#if defined(GPAMELA_UNIX)
C      TSPH(4) = ASIN(SHEL(2)/TSPH(2)) * 180./ACOS(-1.)
       TSPH(4) = ASIN(SHEL(5)/TSPH(2)) * 180./ACOS(-1.)
#endif
#if !defined(GPAMELA_UNIX)
C      TSPH(4) = ASIND(SHEL(2)/TSPH(2))
        TSPH(4) = ASIND(SHEL(5)/TSPH(2))
#endif
*END AB.

*
* 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
*
*AB:
C      SHEL(3) = (ZTRD + TRDB(3) +PAME(3) )/2.
*END: AB.
*
* Add user defined hieght to the external shell
*
*AB:
C      SHEL(3) = SHEL(3) + DZSH
C      SHEI(3) = SHEL(3)
      SHEI(3)=0.
C      ZSHEL = SHEL(3) - PAME(3)
      ZSHEL=ZCAL+CALS(3)+4.+SHEL(1)
C      ZSHEI = SHEI(3) - PAME(3)
      ZSHEI=0.
C      ZTFLA = ZSHEL + SHEL(3) + TFLA(3)
       ZTFLA = ZSHEL + SHEL(1) + TFLA(3)
#if defined(GPAMELA_UNIX)
C      ZTSPH = (2.*SHEI(3)-PAME(3))-TSPH(2)*COS(TSPH(4)/180.*ACOS(-1.))
      ZTSPH=ZSHEL+SHEL(1)-TSPH(2)*COS(TSPH(4)/180.*ACOS(-1.))
#endif
#if !defined(GPAMELA_UNIX)
C      ZTSPH = (2.*SHEI(3)-PAME(3)) - TSPH(2) * COSD(TSPH(4))
    ZTSPH =ZSHEL+SHEL(1)-TSPH(2)*COSD(TSPH(4))
#endif
      ZMSHE=ZSHEL-SHEL(1)-MSHE(3)
      ZBSPH=ZMSHE-MSHE(3)
*END: AB.
C ml: 30/3/06: ZORMG is the offset in z to go from the reference system of
c PAME to the reference system of the tracker (it is the center of the third
c magnetic module). It is obtained so:
c zormg=zcal+cals(3)+zspebp+hormg. zcal (13.05 cm) is the position of the
c center of the calorometer box in PAME, CALS(3) is the half height of the
c calorimeter (10.639 cm), zspebp (2.97 cm) is the distance
c between the bottom of the last alluminum frame of the tracker and the
c bottom of the base plate, HORMG (22.57 cm) is the distance between the
C center of the third magnetic module and the bottom of the tracker
C (22.66 cm is the distance between the center of the third magnetic module
C and the top of the tracker).
      ZORMG=ZCAL+CALS(3)+ZSPEBP+HORMG
c     S11:
c      ZS1=ZPAMS12X+ZORMG+S12X(3)+HMYLTOF+(HPVCMAX-HPVCMIN)/2.
      ZS1=ZCAL+CALS(3)+ZSPEBP+2*SPEB(3)+2*TPLA(3)+HTABLE+DTABLEAL+
     +    S1(3)
      ZTRD=ZCAL+CALS(3)+ZSPEBP+2*SPEB(3)+2*TPLA(3)+TRDB(3)
c now S2 is inside TPLA. So we don't need ZS2, the old position of S2 in PAME:
c In each case, the new position od S2 in PAME has to be calculated as:
C      ZS2=ZCAL+CALS(3)+ZSPEBP+2*SPEB(3)+2*TPCV(3)+S2(3)
C     TPLA (Top Plate):
      ZTPLA=ZCAL+CALS(3)+ZSPEBP+2*SPEB(3)+TPLA(3)
c     tracker:  
      ZSPEC=ZCAL+CALS(3)+ZSPEBP+SPEB(3)
C     CASA:
      ZCAS=ZSPEC-2.1
C     CATA:
c ml: now CATA is in the volume TPLA. The following calculation of ZCAT is
c not necessary in the code. Right to have the position of CATA in PAME: 
c      ZCAT=ZCAL+CALS(3)+ZSPEBP+2*SPEB(3)+2*TPLA(3)-CATA(2)
C     S3 (DTPS3 is the distance between the bottom of the top plate and the
c     bottom of S32):
      ZS3=ZCAL+CALS(3)+DTPS3+HMYLTOF+2*S32X(3)+2.*HMYLTOF+HAIRS3/2.
c     S4:
      ZS4=ZCAL-CALS(3)-S4AL(3)
c     ND:
      ZND=ZCAL-CALS(3)-2*S4AL(3)-NDBO(3)
*
*  Temporary CARD quote by Caf.
*  
      ZCARD = ZTRD
*
* Origin of the spectrometer in the MARS GEANT coordinates system
*
      XYZSOR(1) = 0.
      XYZSOR(2) = 0.
c ml: 30/03/06:
c      XYZSOR(3) = ZSPEC
      XYZSOR(3) = ZORMG
*
* Calculate generation surface & quote based on the S11 dimentions
*
      IF(XYZGEN(1).NE.-1111) THEN
         XGEN = XYZGEN(1)
      ELSE
c         XGEN = -S11(1)
         XGEN=-S1(1)
      ENDIF
      IF(XYZGEN(2).NE.-1111) THEN
         YGEN = XYZGEN(2)
      ELSE
C         YGEN = -S11(2)
         YGEN=-S1(2)
      ENDIF
      IF(XYZGEN(3).NE.-1111) THEN
         ZGEN = XYZGEN(3)
      ELSE
C         ZGEN = ZS11 + 2.
         ZGEN=ZS1+2.
      ENDIF
      IF(XYZGEN(4).NE.-1111) THEN
         XDGEN = XYZGEN(4)
      ELSE
C         XDGEN = 2.*S11(1)
         XDGEN = 2.*S1(1)
      ENDIF
      IF(XYZGEN(5).NE.-1111) THEN
         YDGEN = XYZGEN(5)
      ELSE
C         YDGEN= 2.*S11(2)
         YDGEN= 2.*S1(2)
      ENDIF
      IF(XYZGEN(6).NE.-1111) THEN
         ZDGEN = XYZGEN(6)
      ELSE
         ZDGEN= 0.
      ENDIF
C
C Set max and min generation angles
C
      THMIN = THGEN(1)
      THMAX = THGEN(2)
      COS2G = (THGEN(3).EQ.1.)      
      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.NE.0) 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
C         CALL HCDIR('//FIELD',' ')
C         CALL HRIN(30,9999,0)
C         CALL HPRNT(30)
C         CALL HBNAME(30,' ',0,'$CLEAR')
C         CALL HBNAME(30,'FIELD',F,'$SET')
C         CALL HNOENT(30,IPTOT)
C         IP=0
C         DO IX=1,41
C            DO IY=1,41
C               DO IZ=1,31
C                  IP=IP+1
C                  CALL HGNT(30,IP,IERR)
C                  FX(IX,IY,IZ)=F(1)
C                  FY(IX,IY,IZ)=F(2)
C                  FZ(IX,IY,IZ)=F(3)
C               END DO
C            END DO
C         END DO
C         CALL HREND('FIELD')

C
C fill magnetic field variable needed for interpolation
C
         CALL GPRBI

      ENDIF

C
C READ THE HISTOGRAMS FOR SPE HIT RESOLUTION AND FILL VECTORS
C
      CALL HCDIR('//RESSPE',' ')
      CALL HRIN(0,99999,0)
      DO IK=0,10
         CALL HUNPAK(1000+2*IK,WX(1,IK+1),' ',1)
         CALL HDELET(1000+2*IK)
         CALL HUNPAK(2000+2*IK,WY(1,IK+1),' ',1)
         CALL HDELET(2000+2*IK)
C PREPARE FOR HISTOGRAMS SAMPLING      
         CALL GPRNHPRE(WX(1,IK+1),NDHIRES)
         CALL GPRNHPRE(WY(1,IK+1),NDHIRES)
      ENDDO
      CALL HREND('RESSPE')
C
C prepare the vector for the global strip position to be filled runtime
C
      DO I=1,NPLANES
       DO II=1,NSENSORS
          DO III=1,NSTRIPX
         GLOBSTRIPX(I,II,III)=0.
          ENDDO
          DO III=1,NSTRIPY
          GLOBSTRIPY(I,II,III)=0.
          ENDDO
       ENDDO
      ENDDO 
*
* Set a new mechanism name for Gigantic resonance
*
*      MEC1(IGMEC) = 'GIG '  
*
* Read gigantic resonance cross section
*
      CALL GPGRSIGMA()
*
* 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
