GEANT321/gdraw/gdrack.F

   1 *
   2 * $Id$
   3 *
   4 * $Log$
   5 * Revision 1.1.1.1  1995/10/24 10:20:30  cernlib
   6 * Geant
   7 *
   8 *
   9 #include "geant321/pilot.h"
  10 *CMZ :  3.21/02 29/03/94  15.41.28  by  S.Giani
  11 *-- Author :
  12       SUBROUTINE GDRACK
  13 C.
  14 C.    ******************************************************************
  15 C.    *                                                                *
  16 C.    *       RAY-TRACING                                              *
  17 C.    *       Computation and tracking of all the light rays.          *
  18 C.    *                                                                *
  19 C.    *       Author: S.Giani                                          *
  20 C.    *                                                                *
  21 C.    ******************************************************************
  22 C.
  23 #include "geant321/gcflag.inc"
  24 ********************************************************************************
  25 #include "geant321/gcdraw.inc"
  26 #include "geant321/gconst.inc"
  27 #include "geant321/gcmutr.inc"
  28 #include "geant321/pawc.inc"
  29 ********************************************************************************
  30 #include "geant321/gcbank.inc"
  31 #include "geant321/gccuts.inc"
  32 #include "geant321/gcjloc.inc"
  33 #include "geant321/gckine.inc"
  34 #include "geant321/gcking.inc"
  35 #include "geant321/gcmate.inc"
  36 #include "geant321/gcphys.inc"
  37 #include "geant321/gcparm.inc"
  38 #include "geant321/gcsets.inc"
  39 #include "geant321/gcstak.inc"
  40 #include "geant321/gctmed.inc"
  41 #include "geant321/gctrak.inc"
  42 #include "geant321/gcvolu.inc"
  43 #include "geant321/gcunit.inc"
  44 #include "geant321/gcnum.inc"
  45 #if defined(CERNLIB_USRJMP)
  46 #include "geant321/gcjump.inc"
  47 #endif
  48 ********************************************************************************
  49 #include "geant321/gcfdim.inc"
  50 #include "geant321/gcpixe.inc"
  51 #include "geant321/gcrayt.inc"
  52 #include "geant321/gcvdma.inc"
  53
  54       CHARACTER*4 NAME
  55       DIMENSION VVVMIN(80),VVVMAX(80)
  56       SAVE VVVMIN,VVVMAX,ZIMPRE
  57 ********************************************************************************
  58       COMMON/GCCHAN/LSAMVL
  59       LOGICAL LSAMVL
  60 *
  61       DIMENSION CUTS(10),MECA(5,13)
  62       EQUIVALENCE (CUTS(1),CUTGAM),(MECA(1,1),IPAIR)
  63       SAVE PRECOR
  64 #if !defined(CERNLIB_SINGLE)
  65       PARAMETER (EPSMAC=1.E-6)
  66 #endif
  67 #if defined(CERNLIB_SINGLE)
  68       PARAMETER (EPSMAC=1.E-11)
  69 #endif
  70 C.
  71 C.    ------------------------------------------------------------------
  72          FINVIS=1./ISFILL
  73          LLL=1
  74          ZLN=0.
  75          IF(IMAP.EQ.1)ZIMPRE=LIMPRE
  76          IF(IMAP.EQ.2)THEN
  77            ZRATIO=LIMPRE/ZIMPRE
  78          ENDIF
  79          DO 21 IIIJ=1,LIMPRE
  80            IF(IMAP.EQ.1)THEN
  81             VVVMIN(IIIJ)=10000.
  82             VVVMAX(IIIJ)=0.
  83            ELSEIF(IMAP.EQ.2)THEN
  84             ZLN=ZLN+1
  85             IF((ZLN-ZRATIO).GT.0.001)THEN
  86               ZLN=ZLN-ZRATIO
  87               LLL=LLL+1
  88             ENDIF
  89            ENDIF
  90            VVV=FINVIS
  91            IYYY=1
  92            IF(IIIJ.NE.1)THEN
  93              UUU=UUU+FINVIS
  94              IXXX=IXXX+1
  95            ENDIF
  96           DO 23 JJJI=1,LIMPRE
  97            IF(IMAP.EQ.2)THEN
  98             IF(LLL.EQ.1)THEN
  99              VVMA=VVVMAX(LLL)
 100              VVMI=VVVMIN(LLL)
 101             ELSE
 102              VVMA=MAX(VVVMAX(LLL-1),VVVMAX(LLL))
 103              VVMI=MIN(VVVMIN(LLL-1),VVVMIN(LLL))
 104             ENDIF
 105             IF(((VVV-(VVMA+ZNMAP1)).GT.0.001).OR.
 106      +        ((VVV-(VVMI-ZNMAP1)).LT.-0.001))GOTO 22
 107            ENDIF
 108             IF(IIIJ.EQ.1.AND.JJJI.EQ.1)THEN
 109                MYTRME=NUMED
 110                ime=0
 111 *               print *,vect(1),vect(2),vect(3),'vertex from gtrack'
 112 *               print *,vect(4),vect(5),vect(6),'impulse from gtrack'
 113                GOTO 9
 114             ENDIF
 115                                  XPINTS=ZROTS(1,4)+ZROTS(1,1)*
 116      +                           UUU+ZROTS(1,2)*VVV+ZROTS(1,3)*
 117      +                           ZUV
 118                                  YPINTS=ZROTS(2,4)+ZROTS(2,1)*
 119      +                           UUU+ZROTS(2,2)*VVV+ZROTS(2,3)*
 120      +                           ZUV
 121                                  ZPINTS=ZROTS(3,4)+ZROTS(3,1)*
 122      +                           UUU+ZROTS(3,2)*VVV+ZROTS(3,3)*
 123      +                           ZUV
 124             JON=0
 125             ISSEEN=0
 126             IME=0
 127             SLENG=0.
 128             NUMED=MYTRME
 129             NSTEP=0
 130             INFROM=0
 131             nlevel=1
 132        CALL UCTOH('PERS',IPERS,4,4)
 133        XCOSXS=(SIN(GTHETA*DEGRAD))*(COS(GPHI*DEGRAD))
 134        YCOSYS=(SIN(GTHETA*DEGRAD))*(SIN(GPHI*DEGRAD))
 135        ZCOSZS=COS(GTHETA*DEGRAD)
 136        VDX=XCOSXS
 137        VDY=YCOSYS
 138        VDZ=ZCOSZS
 139        VECT(1)=XPINTS
 140        VECT(2)=YPINTS
 141        VECT(3)=ZPINTS
 142        IF(IPERS.EQ.IPRJ)THEN
 143         CONMOD=1./SQRT(((XPINTS-FPINTX)**2)+((YPINTS-FPINTY)**2)+
 144      +                 ((ZPINTS-FPINTZ)**2))
 145         XCOSXS=-(XPINTS-FPINTX)*CONMOD
 146         YCOSYS=-(YPINTS-FPINTY)*CONMOD
 147         ZCOSZS=-(ZPINTS-FPINTZ)*CONMOD
 148        ENDIF
 149        VECT(4) = -XCOSXS
 150        VECT(5) = -YCOSYS
 151        VECT(6) = -ZCOSZS
 152        IF(VECT(1).LE.XCUT.OR.VECT(2).LE.YCUT.OR.VECT(3).LE.ZCUT)THEN
 153          CALL GTMEDI (VECT, NUMED)
 154        ENDIF
 155  9    CONTINUE
 156
 157       ISTOP = 0
 158       EPSCUR = EPSMAC
 159       NSTOUT = 0
 160       INWOLD = 0
 161       LSAMVL = .FALSE.
 162       NUMOLD=0
 163 *
 164 * *** Check validity of tracking medium and material parameters
 165 *
 166    10 IF (NUMED.NE.NUMOLD) THEN
 167          NUMOLD = NUMED
 168          IUPD   = 0
 169          JTM    = LQ(JTMED- NUMED)
 170          EPSIL    = Q(JTM + 13)
 171          PRECOR   = MIN(0.1*EPSIL, 0.0010)
 172 *
 173          NMAT     = Q(JTM + 6)
 174          JMA   = LQ(JMATE-NMAT)
 175          DENS = Q(JMA +8)
 176       ENDIF
 177         IF(ISCOLO.EQ.-10.OR.IMYSE.EQ.0)ISCOLO=MOD(NUMED,6)+2
 178         IF(ISLSTY.EQ.-10.OR.IMYSE.EQ.0)ISLSTY=4
 179         IF(ISSEEN.EQ.-10.OR.IMYSE.EQ.0)THEN
 180          IF(IME.EQ.1)THEN
 181          IF(DENS.LT.0.00130)THEN
 182           ISSEEN=0
 183          ELSE
 184           ISSEEN=1
 185          ENDIF
 186          ENDIF
 187         ENDIF
 188         IF(IME.EQ.0)IME=1
 189 *
 190       IF(LSAMVL) THEN
 191 *
 192 *       If now the particle is entering in the same volume where
 193 *       it was exiting from last step, and if it has done this for
 194 *       more than 5 times, we decrease the precision of tracking
 195          NSTOUT=NSTOUT+1
 196          IF(MOD(NSTOUT,5).EQ.0) THEN
 197             EPSCUR=NSTOUT*EPSMAC
 198 *            WRITE(CHMAIL,10000)ITRA,ISTAK,NTMULT,NAPART
 199 *10000          FORMAT(' *** GTRACK *** Boundary loop: track ',
 200 *     +         I4,' stack ',I4,' NTMULT ',I5,1X,5A4)
 201 *            CALL GMAIL(1,0)
 202 *            WRITE(CHMAIL,10100) EPSCUR
 203 *10100          FORMAT('                Precision now set to ',G10.3)
 204 *            CALL GMAIL(0,1)
 205          ENDIF
 206       ELSE
 207          NSTOUT = 0
 208          EPSCUR = EPSMAC
 209       ENDIF
 210 *
 211       INWVOL = 1
 212 *
 213 * *** Compute SET and DET number if volume is sensitive
 214 *
 215       IF (JSET.GT.0) CALL GFINDS
 216 *
 217 *    Clear step dependent variables
 218 *
 219    80 NMEC   = 0
 220       STEP   = 0.
 221       DESTEL = 0.
 222       DESTEP = 0.
 223       NGKINE = 0
 224       IGNEXT = 0
 225       INWOLD = INWVOL
 226       PREC   = MAX(PRECOR,MAX(ABS(VECT(1)),ABS(VECT(2)),
 227      +                        ABS(VECT(3)),SLENG)*EPSCUR)
 228 *
 229 *     Give control to user at entrance of volume (INWVOL=1)
 230 *
 231       IF (INWVOL.EQ.1) THEN
 232          CALL GDSTEP
 233          IF (ISTOP.NE.0) GO TO 22
 234          INWVOL = 0
 235       ENDIF
 236 *
 237 * *** Propagate particle up to next volume boundary or end of track
 238 *
 239       INGOTO = 0
 240       NLEVIN = NLEVEL
 241       IF (IPARAM.NE.0) THEN
 242          IF (GEKIN.LE.PACUTS(ITRTYP)) THEN
 243             NMEC = NMEC+1
 244             LMEC(NMEC) = 26
 245             ISTOP = 2
 246 #if !defined(CERNLIB_USRJMP)
 247             CALL GUPARA
 248 #endif
 249 #if defined(CERNLIB_USRJMP)
 250             CALL JUMPT0(JUPARA)
 251 #endif
 252             GO TO 90
 253          ENDIF
 254       ENDIF
 255
 256          CALL GDNINO
 257
 258       STLOSS=STEP
 259 *
 260 *     Check for possible endless loop
 261 *
 262    90 NSTEP = NSTEP +1
 263
 264       IF (NSTEP.GT.MAXNST) THEN
 265          IF (ISTOP.EQ.0) THEN
 266             ISTOP = 99
 267             NMEC  = NMEC +1
 268             LMEC(NMEC) = 30
 269 *            WRITE(CHMAIL,10200) MAXNST
 270 *            CALL GMAIL(1,0)
 271 *            WRITE(CHMAIL,10300)ITRA,ISTAK,NTMULT,(NAPART(I),I=1,5),
 272 *     +      TOFG*1.E9
 273 *            CALL GMAIL(0,1)
 274 *10200       FORMAT(' *** GTRACK *** More than ',I6,
 275 *     +             ' steps, tracking abandoned!')
 276 *10300       FORMAT('                Track ',I4,' stack ',I4,' NTMULT ',
 277 *     +             I5,1X,5A4,1X,'Time of flight ',F10.3,' ns')
 278          ENDIF
 279       ENDIF
 280 *
 281 * *** Give control to user at end of each tracking step
 282 *
 283       SAFETY = SAFETY -STEP
 284       CALL GDSTEP
 285 *
 286       IF (ISTOP.NE.0) GO TO 22
 287 *
 288 *      Renormalize direction cosines
 289 *
 290       CMOD = 1./SQRT(VECT(4)*VECT(4)+VECT(5)*VECT(5)+VECT(6)*VECT(6))
 291       VECT(4) = VECT(4)*CMOD
 292       VECT(5) = VECT(5)*CMOD
 293       VECT(6) = VECT(6)*CMOD
 294 *
 295       IF (INWVOL.EQ.0) GO TO 80
 296 *
 297       IF (NJTMAX.GT.0) THEN
 298          CALL GSTRAC
 299          IF (NLEVIN.EQ.0) GO TO 100
 300          GO TO 22
 301       ELSE
 302          IF (NLEVIN.GE.NLEVEL) THEN
 303             INFROM = 0
 304          ELSE
 305             IF (NLEVIN.EQ.0) GO TO 100
 306             INFROM = LINDEX(NLEVIN+1)
 307          ENDIF
 308          IF (NLEVIN.NE.NLEVEL) INGOTO = 0
 309          NLEVEL = NLEVIN
 310 *
 311          CALL GTMEDI (VECT, NUMED)
 312          IF (NUMED.NE.0) THEN
 313             SAFETY = 0.
 314             CALL UHTOC(NAMES(NLEVEL),4,NAME,4)
 315 *            print *,NAME
 316            IF(IMYSE.EQ.1)THEN
 317             CALL GFIND(NAME,'SEEN',ISSEEN)
 318             CALL GFIND(NAME,'COLO',ISCOLO)
 319             CALL GFIND(NAME,'LSTY',ISLSTY)
 320            ENDIF
 321 *            print *,isseen,iscolo,islsty
 322             GO TO 10
 323          ELSE
 324             ISSEEN=1
 325          ENDIF
 326       ENDIF
 327 *
 328 *     Track outside setup, give control to user (INWVOL=3)
 329 *
 330   100 INWVOL = 3
 331       ISTOP  = 1
 332       ISET   = 0
 333       IDET   = 0
 334       NMEC   = 0
 335       STEP   = 0.
 336       DESTEL = 0.
 337       DESTEP = 0.
 338       NGKINE = 0
 339       NLCUR  = NLEVEL
 340       NLEVEL = 1
 341       CALL GDSTEP
 342       NLEVEL = NLCUR
 343  22        CONTINUE
 344           IF(IMAP.EQ.1)THEN
 345            IF(JON.EQ.1)THEN
 346              IF(VVV.LT.VVVMIN(IIIJ))THEN
 347                VVVMIN(IIIJ)=VVV
 348              ELSEIF(VVV.GT.VVVMAX(IIIJ))THEN
 349                VVVMAX(IIIJ)=VVV
 350              ENDIF
 351            ENDIF
 352           ENDIF
 353            VVV=VVV+FINVIS
 354            IYYY=IYYY+1
 355    23     CONTINUE
 356    21    CONTINUE
 357 *                                                             END GTRACK
 358   999 END