GEANT321/gtrak/gtmedi.F

   1 *
   2 * $Id$
   3 *
   4 * $Log$
   5 * Revision 1.1.1.1  1995/10/24 10:21:44  cernlib
   6 * Geant
   7 *
   8 *
   9 #include "geant321/pilot.h"
  10 #if !defined(CERNLIB_OLD)
  11 *CMZ :  3.21/02 29/03/94  15.22.08  by  S.Giani
  12 *-- Author :
  13       SUBROUTINE GTMEDI (X, NUMED)
  14 C.
  15 C.    ******************************************************************
  16 C.    *                                                                *
  17 C.    *   Finds in which volume/medium the point X is, and updates the *
  18 C.    *    common /GCVOLU/ and the structure JGPAR accordingly.        *
  19 C.    *                                                                *
  20 C.    *   NUMED returns the tracking medium number, or 0 if point is   *
  21 C.    *         outside the experimental setup.                        *
  22 C.    *                                                                *
  23 C.    *   Note : For INWVOL = 2, INFROM set to a positive number is    *
  24 C.    *      interpreted by GTMEDI as the number IN of the content     *
  25 C.    *      just left by the current track within the mother volume   *
  26 C.    *      where the point X is assumed to be.                       *
  27 C.    *                                                                *
  28 C.    *   Note : INFROM is set correctly by this routine but it is     *
  29 C.    *      used on entrance only in the case GSNEXT has been called  *
  30 C.    *      by the user. In other words the value of INFROM received  *
  31 C.    *      on entrance is not considered necessarily valid. This     *
  32 C.    *      assumption has been made for safety. A wrong value of     *
  33 C.    *      INFROM can cause wrong tracking.                          *
  34 C.    *                                                                *
  35 C.    *   Called by : GTRACK                                           *
  36 C.    *   Authors   : S.Banerjee, R.Brun, F.Bruyant, A.McPherson       *
  37 C.    *               S.Giani.                                         *
  38 C.    *                                                                *
  39 C.    *   Modified by S.Giani (1993) to perform the search according   *
  40 C.    *    to the new 'virtual divisions' algorithm and to build the   *
  41 C.    *    stack of the 'possible overlapping volumes' in the case of  *
  42 C.    *    MANY volumes. Any kind of boolean operation is now possible.*
  43 C.    *    Divisions along arbitrary axis are now possible.            *
  44 C.    *                                                                *
  45 C.    ******************************************************************
  46 C.
  47 #include "geant321/gcflag.inc"
  48 #include "geant321/gckine.inc"
  49 #include "geant321/gcbank.inc"
  50 #include "geant321/gcvolu.inc"
  51 #include "geant321/gctrak.inc"
  52 #if defined(CERNLIB_USRJMP)
  53 #include "geant321/gcjump.inc"
  54 #endif
  55 #include "geant321/gchvir.inc"
  56 #include "geant321/gcvdma.inc"
  57       COMMON/GCCHAN/LSAMVL
  58       LOGICAL LSAMVL
  59 C.
  60       CHARACTER*4 NAME
  61       DIMENSION  X(*)
  62       REAL       XC(6), XT(3)
  63       LOGICAL    BTEST
  64 C.
  65 C.    ------------------------------------------------------------------
  66 *
  67       nvmany=0
  68       nfmany=0
  69       new2fl=0
  70       neufla=0
  71       if(raytra.eq.1.)then
  72         JVO  = LQ(JVOLUM-LVOLUM(NLEVEL))
  73         if(ingoto.eq.-1.and.q(jvo+3).lt.0)then
  74           neufla=1
  75         elseif(ingoto.eq.0)then
  76               call ggperp(x,veccos,ierr)
  77               veccos(1)=-veccos(1)
  78               veccos(2)=-veccos(2)
  79               veccos(3)=-veccos(3)
  80               if(ierr.eq.1)then
  81                 veccos(1)=1.
  82                 veccos(2)=0.
  83                 veccos(3)=0.
  84               endif
  85         endif
  86       endif
  87 *
  88 * SECTION I: The /GCVOLU/ table contains the initial guess for a path
  89 *            in the geometry tree on which X may be found.  Look along this
  90 *            path until X is found inside.  This is the starting position.
  91 *            If this is an ONLY volume with no daughters, we are done;
  92 *            otherwise reset search record variables, proceed to section II.
  93 *
  94 * *** Check if point is in current volume
  95 *
  96       INFR = 0
  97       INGT = 0
  98       JVIN = 0
  99 *
 100 * *** LSAMVL is a logical variable that indicates whether we are still
 101 * *** in the current volume or not. It is used in GTRACK to detect
 102 * *** precision problems.
 103       LSAMVL = .TRUE.
 104 C*****  Code Expanded From Routine:  GTRNSF
 105 C
 106   100 IF (GRMAT(10,NLEVEL) .EQ. 0.) THEN
 107          XC(1) = X(1) - GTRAN(1,NLEVEL)
 108          XC(2) = X(2) - GTRAN(2,NLEVEL)
 109          XC(3) = X(3) - GTRAN(3,NLEVEL)
 110 *
 111       ELSE
 112          XL1 = X(1) - GTRAN(1,NLEVEL)
 113          XL2 = X(2) - GTRAN(2,NLEVEL)
 114          XL3 = X(3) - GTRAN(3,NLEVEL)
 115          XC(1) = XL1*GRMAT(1,NLEVEL) + XL2*GRMAT(2,NLEVEL) + XL3*
 116      +      GRMAT(3,NLEVEL)
 117          XC(2) = XL1*GRMAT(4,NLEVEL) + XL2*GRMAT(5,NLEVEL) + XL3*
 118      +      GRMAT(6,NLEVEL)
 119          XC(3) = XL1*GRMAT(7,NLEVEL) + XL2*GRMAT(8,NLEVEL) + XL3*
 120      +      GRMAT(9,NLEVEL)
 121
 122       ENDIF
 123       xc(4)=0.
 124       xc(5)=0.
 125       xc(6)=0.
 126 C*****  End of Code Expanded From Routine:  GTRNSF
 127 *
 128       JVO  = LQ(JVOLUM-LVOLUM(NLEVEL))
 129 *
 130 * Note: At entry the variable INGOTO may contain the index of a volume
 131 * contained within the current one at NLEVEL.  If so, begin by checking
 132 * if X lies inside.  This improves the search speed over that of GMEDIA.
 133 *
 134       NIN = Q(JVO+3)
 135       if(raytra.eq.1..and.imyse.eq.1)then
 136             CALL UHTOC(NAMES(NLEVEL),4,NAME,4)
 137             CALL GFIND(NAME,'SEEN',ISSEEN)
 138             if(isseen.eq.-2.or.isseen.eq.-1)goto 189
 139       endif
 140       IF ((INGOTO.LE.0).OR.(INGOTO.GT.NIN)) THEN
 141          INGOTO = 0
 142       ELSE
 143 *
 144 * ***   Entrance in content INGOTO predicted by GTNEXT
 145 *
 146          JIN  = LQ(JVO-INGOTO)
 147          IVOT = Q(JIN+2)
 148          JVOT = LQ(JVOLUM-IVOT)
 149          JPAR = LQ(JGPAR-NLEVEL-1)
 150 *
 151          IROTT = Q(JIN+4)
 152 C*****  Code Expanded From Routine:  GITRAN
 153 C.
 154 C.    ------------------------------------------------------------------
 155 C.
 156          IF (IROTT .EQ. 0) THEN
 157             XT(1) = XC(1) - Q(5+JIN)
 158             XT(2) = XC(2) - Q(6+JIN)
 159             XT(3) = XC(3) - Q(7+JIN)
 160 *
 161          ELSE
 162             XL1 = XC(1) - Q(5+JIN)
 163             XL2 = XC(2) - Q(6+JIN)
 164             XL3 = XC(3) - Q(7+JIN)
 165             JR = LQ(JROTM-IROTT)
 166             XT(1) = XL1*Q(JR+1) + XL2*Q(JR+2) + XL3*Q(JR+3)
 167             XT(2) = XL1*Q(JR+4) + XL2*Q(JR+5) + XL3*Q(JR+6)
 168             XT(3) = XL1*Q(JR+7) + XL2*Q(JR+8) + XL3*Q(JR+9)
 169 *
 170          ENDIF
 171 C*****  End of Code Expanded From Routine:  GITRAN
 172 *
 173 *   *   Check if point is in content
 174 *
 175          CALL GINME (XT, Q(JVOT+2), Q(JPAR+1), IYES)
 176          IF (IYES.NE.0) THEN
 177 *
 178 *          If so, prepare information for volume retrieval, and return
 179 *
 180             LSAMVL = .FALSE.
 181             NL1 = NLEVEL +1
 182             LVOLUM(NL1) = IVOT
 183             NAMES(NL1)  = IQ(JVOLUM+IVOT)
 184             NUMBER(NL1) = Q(JIN+3)
 185             LINDEX(NL1) = INGOTO
 186             LINMX(NL1)  = Q(JVO+3)
 187             GONLY(NL1)  = Q(JIN+8)
 188             IF (LQ(LQ(JVOLUM-IVOT)).EQ.0) THEN
 189                NLDEV(NL1) = NLDEV(NLEVEL)
 190             ELSE
 191                NLDEV(NL1) = NL1
 192             ENDIF
 193             CALL GTRMUL (GTRAN(1,NLEVEL), GRMAT(1,NLEVEL), Q(JIN+5),
 194      +                   IROTT, GTRAN(1,NL1), GRMAT(1,NL1))
 195             NLEVEL = NL1
 196             XC(1) = XT(1)
 197             XC(2) = XT(2)
 198             XC(3) = XT(3)
 199             JVO = JVOT
 200             INFROM = 0
 201             if(raytra.eq.1.)then
 202               call ggperp(x,veccos,ierr)
 203               if(ierr.eq.1)then
 204                 veccos(1)=1.
 205                 veccos(2)=0.
 206                 veccos(3)=0.
 207               endif
 208             endif
 209             GO TO 190
 210          ENDIF
 211       ENDIF
 212 *
 213 * End of INGOTO processing
 214 *
 215  189  JPAR = LQ(JGPAR-NLEVEL)
 216       CALL GINME (XC, Q(JVO+2), Q(JPAR+1), IYES)
 217       IF (IYES.EQ.0) THEN
 218 *
 219 *  **   Point not in current volume, go up the tree
 220 *
 221          LSAMVL = .FALSE.
 222          INGOTO = 0
 223          IF (NLEVEL.GT.1) THEN
 224             NLEVEL = NLEVEL -1
 225             JVO = LQ(JVOLUM-LVOLUM(NLEVEL))
 226             NIN = Q(JVO+3)
 227             IF(NIN.GT.0) THEN
 228                INFROM=LINDEX(NLEVEL+1)
 229             ELSE
 230                INFROM=0
 231             ENDIF
 232             INFR = INFROM
 233             GO TO 100
 234          ELSE
 235 *
 236 *   *      Point is outside setup
 237 *
 238             NUMED = 0
 239             GO TO 999
 240          ENDIF
 241       ELSE
 242 *
 243 *   *      Point in current volume but not in INGOTO. We block the
 244 *   *      corresponding volume
 245 *
 246          IF (INGOTO.GT.0) THEN
 247             INGT = INGOTO
 248             JIN = LQ(JVO-INGOTO)
 249             IQ(JIN) = IBSET(IQ(JIN),4)
 250          ENDIF
 251       ENDIF
 252 *
 253 *   *      Found a volume up the tree which contains our point. We block
 254 *   *      the branch we came up from.
 255 *
 256       IF(INFR.GT.0) THEN
 257          JIN=LQ(JVO-INFR)
 258          IQ(JIN) = IBSET(IQ(JIN),4)
 259          JVIN = JIN
 260       ENDIF
 261 *
 262 *  **   Point is in current volume
 263 *
 264   190 INGOTO = 0
 265       NLMIN = NLEVEL
 266       IF (INWVOL.NE.2) INFROM = 0
 267       NLMANY = 0
 268 *
 269 * SECTION II: X is found inside current node at NLEVEL in /GCVOLU/.
 270 *             Search all contents recursively for any containing X.
 271 *             Take the first one found, if any, and continue at that
 272 *             level, incrementing NLEVEL and extending /GCVOLU/ tables.
 273 *             This is continued until a level is reached where X is not
 274 *             found in any of the contents, or there are no contents.
 275 * Note: Since Section II is re-entered from Section III, a blocking word
 276 * is used to mark those contents already checked.  Upon exit from Section
 277 * II, these blocking words are cleared at NLEVEL, but may remain set in
 278 * levels between NLEVEL-1 and NLMIN, if any.  They must be cleared at exit.
 279 *
 280 *  **  Check contents, if any
 281 *
 282   200 JVO = LQ(JVOLUM-LVOLUM(NLEVEL))
 283       NIN = Q(JVO+3)
 284       if(raytra.eq.1..and.imyse.eq.1)then
 285             CALL UHTOC(NAMES(NLEVEL),4,NAME,4)
 286             CALL GFIND(NAME,'SEEN',ISSEEN)
 287             if(isseen.eq.-2.or.isseen.eq.-1)goto 300
 288       endif
 289 *
 290 *   *   Case with no contents
 291 *
 292       IF (NIN.EQ.0) THEN
 293          GO TO 300
 294 *
 295 *   *   Case with contents defined by division
 296 *
 297       ELSEIF (NIN.LT.0) THEN
 298          CALL GMEDIV (JVO, IN, XC, 1)
 299          IF (IN.GT.0) THEN
 300             if(raytra.eq.1..and.neufla.eq.1)then
 301               neufla=0
 302               call ggperp(x,veccos,ierr)
 303               if(ierr.eq.1)then
 304                 veccos(1)=1.
 305                 veccos(2)=0.
 306                 veccos(3)=0.
 307               endif
 308             endif
 309             INFROM = 0
 310             INFR   = 0
 311             INGT   = 0
 312             LSAMVL = .FALSE.
 313             GO TO 200
 314          ENDIF
 315 *
 316 *   *  Case with contents positioned
 317 *
 318       ELSE
 319        if(nin.gt.1)then
 320         clmoth=q(jvirt+4*(LVOLUM(NLEVEL)-1)+3)
 321         chmoth=q(jvirt+4*(LVOLUM(NLEVEL)-1)+4)
 322         ndivto=q(jvirt+4*(LVOLUM(NLEVEL)-1)+2)
 323         iaxis =q(jvirt+4*(LVOLUM(NLEVEL)-1)+1)
 324         if(iaxis.le.3)then
 325           ivdiv=((xc(iaxis)-clmoth)*ndivto/(chmoth-clmoth))+1
 326             if(ivdiv.lt.1)then
 327               ivdiv=1
 328             elseif(ivdiv.gt.ndivto)then
 329               ivdiv=ndivto
 330             endif
 331         else
 332           call gfcoor(xc,iaxis,cx)
 333           if(iaxis.eq.6)then
 334             if((cx-clmoth).lt.-1.)then
 335               cx=cx+360.
 336             elseif((cx-chmoth).gt.1.)then
 337               cx=cx-360.
 338             endif
 339             if(cx.gt.chmoth)then
 340               cx=chmoth
 341             elseif(cx.lt.clmoth)then
 342               cx=clmoth
 343             endif
 344           endif
 345           ivdiv=((cx-clmoth)*ndivto/(chmoth-clmoth))+1
 346             if(ivdiv.lt.1)then
 347               ivdiv=1
 348             elseif(ivdiv.gt.ndivto)then
 349               ivdiv=ndivto
 350             endif
 351         endif
 352         jvdiv=lq(jvirt-LVOLUM(NLEVEL))
 353         iofset=iq(jvdiv+ivdiv)
 354         ncont=iq(jvdiv+iofset+1)
 355         jcont=jvdiv+iofset+1
 356         if(ncont.eq.0)goto 260
 357        else
 358          JCONT = LQ(JVO-NIN-1)+1
 359          NCONT = 1
 360        endif
 361 *
 362 *     For each selected content in turn, check if point is inside
 363 *
 364          DO 259 ICONT=1,NCONT
 365            if(nin.eq.1)then
 366             in=1
 367            else
 368             IN = IQ(JCONT+ICONT)
 369            endif
 370             IF(IN.EQ.0) THEN
 371 *
 372 *     If the value IQ(JCONT+ICONT)=0 then we are back in the mother.
 373 *     So jump to 260, the search is finished. Clean-up should be done
 374 *     only up to ICONT-1, so we set:
 375 *
 376                NCONT=ICONT-1
 377                GOTO 260
 378             ELSE
 379             JIN = LQ(JVO-IN)
 380             IF (.NOT.BTEST(IQ(JIN),4)) THEN
 381                CALL GMEPOS (JVO, IN, XC, 1)
 382                IF (IN.GT.0) THEN
 383                   new2fl=0
 384                   IF (GONLY(NLEVEL).NE.0.) THEN
 385                     NLMANY = 0
 386                     nvmany = 0
 387                     nfmany = 0
 388                   ENDIF
 389                   INFROM = 0
 390                   INGT   = 0
 391                   INFR   = 0
 392                   LSAMVL = .FALSE.
 393                   GO TO 200
 394                ELSE
 395                   IQ(JIN) = IBSET(IQ(JIN),4)
 396                ENDIF
 397             ENDIF
 398             ENDIF
 399   259    CONTINUE
 400 *
 401   260    IF(NCONT.EQ.NIN) THEN
 402          DO 268 IN=1,NIN
 403             JIN = LQ(JVO-IN)
 404             IQ(JIN) = IBCLR(IQ(JIN),4)
 405   268    CONTINUE
 406          ELSE
 407          DO 269 ICONT=1,NCONT
 408            if(nin.eq.1)then
 409             in=1
 410            else
 411             IN  = IQ(JCONT+ICONT)
 412            endif
 413             JIN = LQ(JVO-IN)
 414             IQ(JIN) = IBCLR(IQ(JIN),4)
 415   269    CONTINUE
 416          IF(INFR.NE.0) THEN
 417             JIN = LQ(JVO-INFR)
 418             IQ(JIN) = IBCLR(IQ(JIN),4)
 419             INFR = 0
 420          ENDIF
 421          IF(INGT.NE.0) THEN
 422             JIN = LQ(JVO-INGT)
 423             IQ(JIN) = IBCLR(IQ(JIN),4)
 424             INGT = 0
 425          ENDIF
 426          ENDIF
 427 *
 428       ENDIF
 429 *
 430 * SECTION III: X is found at current node (NLEVEL in /GCVOLU) but not in
 431 *              any of its contents, if any.  If this is a MANY volume,
 432 *              save it as a candidate best-choice, and continue the search
 433 *              by backing up the tree one node and proceed to Section II.
 434 *              If this is an ONLY volume, proceed to Section IV.
 435 *
 436 * *** Point is in current volume/medium, and not in any content
 437 *
 438   300 IF (GONLY(NLEVEL).EQ.0.) THEN
 439 *
 440 *  **   Lowest level is 'NOT ONLY'
 441 *
 442          IF (NLEVEL.GT.NLMANY) THEN
 443             CALL GSCVOL
 444             NLMANY = NLEVEL
 445             nfmany=nvmany+1
 446          ENDIF
 447          if(new2fl.eq.0)then
 448             nvmany=nvmany+1
 449             manyle(nvmany)=nlevel
 450             do 401 i = 1,nlevel
 451               manyna(nvmany,i)=names(i)
 452               manynu(nvmany,i)=number(i)
 453  401        continue
 454          endif
 455 *
 456 *   *   Go up the tree up to a volume with positioned contents
 457 *
 458          new2fl=-1
 459   310    INFROM = LINDEX(NLEVEL)
 460          NLEVEL = NLEVEL -1
 461          JVO    = LQ(JVOLUM-LVOLUM(NLEVEL))
 462          NIN    = Q(JVO+3)
 463          IF (NIN.LT.0) GO TO 310
 464 *
 465 C*****  Code Expanded From Routine:  GTRNSF
 466 C
 467          IF (GRMAT(10,NLEVEL) .EQ. 0.) THEN
 468             XC(1) = X(1) - GTRAN(1,NLEVEL)
 469             XC(2) = X(2) - GTRAN(2,NLEVEL)
 470             XC(3) = X(3) - GTRAN(3,NLEVEL)
 471 *
 472          ELSE
 473             XL1 = X(1) - GTRAN(1,NLEVEL)
 474             XL2 = X(2) - GTRAN(2,NLEVEL)
 475             XL3 = X(3) - GTRAN(3,NLEVEL)
 476             XC(1) = XL1*GRMAT(1,NLEVEL) + XL2*GRMAT(2,NLEVEL) +
 477      +      XL3* GRMAT(3,NLEVEL)
 478             XC(2) = XL1*GRMAT(4,NLEVEL) + XL2*GRMAT(5,NLEVEL) +
 479      +      XL3* GRMAT(6,NLEVEL)
 480             XC(3) = XL1*GRMAT(7,NLEVEL) + XL2*GRMAT(8,NLEVEL) +
 481      +      XL3* GRMAT(9,NLEVEL)
 482 *
 483          ENDIF
 484 C*****  End of Code Expanded From Routine:  GTRNSF
 485 *
 486          INFR = INFROM
 487          JIN = LQ(JVO-INFROM)
 488          IQ(JIN) = IBSET(IQ(JIN),4)
 489          NLMIN = MIN(NLEVEL,NLMIN)
 490          GO TO 200
 491       ENDIF
 492 *
 493 * SECTION IV: This is the end of the search.  The current node (NLEVEL
 494 *             in /GCVOLU/) is the lowest ONLY volume in which X is found.
 495 *             If X was also found in any of its contents, they are MANY
 496 *             volumes: the best-choice is the one among them at the greatest
 497 *             level in the tree, and it is stored.  Otherwise the current
 498 *             volume is the solution.  Before exit, all of the blocking
 499 *             words leftover in the tree must be reset to zero.
 500 * Note: A valid structure is assumed, in which no ONLY volumes overlap.
 501 * If this rule is violated, or if a daughter is not entirely contained
 502 * within the mother volume, the results are unpredictable.
 503 *
 504       DO 419 NL=NLMIN,NLEVEL-1
 505          JVO = LQ(JVOLUM-LVOLUM(NL))
 506          NIN = Q(JVO+3)
 507          DO 418 IN=1,NIN
 508             JIN = LQ(JVO-IN)
 509             IQ(JIN) = IBCLR(IQ(JIN),4)
 510   418    CONTINUE
 511   419 CONTINUE
 512 *
 513       if(nlmany.eq.0)then
 514         nvmany=0
 515         nfmany=0
 516       endif
 517       IF (NLMANY.GT.0) CALL GFCVOL
 518       JVO = LQ(JVOLUM-LVOLUM(NLEVEL))
 519       IF(JVIN.NE.0) IQ(JVIN) = IBCLR(IQ(JVIN),4)
 520       NUMED = Q(JVO+4)
 521 *                                                             END GTMEDI
 522   999 IF(JGSTAT.NE.0) CALL GFSTAT(4)
 523       END
 524 #endif