GEANT321/ggeom/gmedia.F

   1 *
   2 * $Id$
   3 *
   4 * $Log$
   5 * Revision 1.3  1998/02/09 16:48:33  japost
   6 *   Simone's fix for MANY volumes in gmedia.
   7 *
   8 * Revision 1.2  1996/09/30 14:25:07  ravndal
   9 * Windows NT related modifications
  10 *
  11 * Revision 1.1.1.1  1995/10/24 10:20:51  cernlib
  12 * Geant
  13 *
  14 *
  15 #include "geant321/pilot.h"
  16 #if !defined(CERNLIB_OLD)
  17 *CMZ :  3.21/02 29/03/94  15.24.17  by  S.Giani
  18 *-- Author :
  19       SUBROUTINE GMEDIA (X, NUMED)
  20 C.
  21 C.    ******************************************************************
  22 C.    *                                                                *
  23 C.    *   Finds in which volume/medium the point X is, and updates the *
  24 C.    *    common /GCVOLU/ and the structure JGPAR accordingly.        *
  25 C.    *                                                                *
  26 C.    *   NUMED returns the tracking medium number, or 0 if point is   *
  27 C.    *         outside the experimental setup.                        *
  28 C.    *                                                                *
  29 C.    *   Called by :  GTREVE, GLTRAC, 'User'                          *
  30 C.    *   Authors   : R.Brun, F.Bruyant, A.McPherson                   *
  31 C.    *               S.Giani.                                         *
  32 C.    *                                                                *
  33 C.    *   Modified by S.Giani (1993) to perform the search according   *
  34 C.    *    to the new 'virtual divisions' algorithm and to build the   *
  35 C.    *    stack of the 'possible overlapping volumes' in the case of  *
  36 C.    *    MANY volumes. Any kind of boolean operation is now possible.*
  37 C.    *    Divisions along arbitrary axis are now possible.            *
  38 C.    *                                                                *
  39 C.    ******************************************************************
  40 C.
  41 #include "geant321/gcflag.inc"
  42 #include "geant321/gckine.inc"
  43 #include "geant321/gcbank.inc"
  44 #include "geant321/gcvolu.inc"
  45 #include "geant321/gctrak.inc"
  46 #if defined(CERNLIB_USRJMP)
  47 #include "geant321/gcjump.inc"
  48 #endif
  49 #include "geant321/gcvdma.inc"
  50 #include "geant321/gchvir.inc"
  51 C.
  52       DIMENSION  X(*)
  53       REAL       XC(6)
  54       LOGICAL    BTEST
  55       CHARACTER*4 NAME
  56 C.
  57 C.    ------------------------------------------------------------------
  58 *
  59       nvmany=0
  60       nfmany=0
  61       new2fl=0
  62 *
  63       IF (NLEVEL.EQ.0) CALL GMEDIN
  64 *
  65 * SECTION I: The /GCVOLU/ table contains the initial guess for a path
  66 *            in the geometry tree on which X may be found.  Look along this
  67 *            path until X is found inside.  This is the starting position.
  68 *            If this is an ONLY volume with no daughters, we are done;
  69 *            otherwise reset search record variables, proceed to section II.
  70 *
  71 *            The information contained in INFROM has to be invalidated
  72 *            because it has no meaning for the subsequent tracking. INFR
  73 *            is a local variable used to optimise the search in the
  74 *            geometry tree.
  75 *
  76       INFROM = 0
  77 *
  78 * *** Check if point is in current volume
  79 *
  80       INFR   = 0
  81       JVIN   = 0
  82 C*****  Code Expanded From Routine:  GTRNSF
  83 C
  84   100 IF (GRMAT(10,NLEVEL) .EQ. 0.) THEN
  85          XC(1) = X(1) - GTRAN(1,NLEVEL)
  86          XC(2) = X(2) - GTRAN(2,NLEVEL)
  87          XC(3) = X(3) - GTRAN(3,NLEVEL)
  88 *
  89       ELSE
  90          XL1 = X(1) - GTRAN(1,NLEVEL)
  91          XL2 = X(2) - GTRAN(2,NLEVEL)
  92          XL3 = X(3) - GTRAN(3,NLEVEL)
  93          XC(1) = XL1*GRMAT(1,NLEVEL) + XL2*GRMAT(2,NLEVEL) + XL3*GRMAT(3
  94      +      ,NLEVEL)
  95          XC(2) = XL1*GRMAT(4,NLEVEL) + XL2*GRMAT(5,NLEVEL) + XL3*GRMAT(6
  96      +      ,NLEVEL)
  97          XC(3) = XL1*GRMAT(7,NLEVEL) + XL2*GRMAT(8,NLEVEL) + XL3*GRMAT(9
  98      +      ,NLEVEL)
  99
 100       ENDIF
 101       xc(4)=0.
 102       xc(5)=0.
 103       xc(6)=0.
 104 C*****  End of Code Expanded From Routine:  GTRNSF
 105 *
 106       JVO  = LQ(JVOLUM-LVOLUM(NLEVEL))
 107       JPAR = LQ(JGPAR-NLEVEL)
 108       CALL GINME (XC, Q(JVO+2), Q(JPAR+1), IYES)
 109       IF (IYES.EQ.0) THEN
 110 *
 111 *  **   Point not in current volume, go up the tree
 112 *
 113          IF (NLEVEL.GT.1) THEN
 114             NLEVEL = NLEVEL -1
 115             JVO  = LQ(JVOLUM-LVOLUM(NLEVEL))
 116             NIN = Q(JVO+3)
 117             IF(NIN.GT.0) THEN
 118 *
 119 *       Do not set INFR whne going up the tree. GMEDIA can be called
 120 *       by the user and it  should not assume  that the previous
 121 *       position has something to do with the current search. INFR
 122 *       is otherwise useful when searching in a 'MANY' volume
 123 *       configuration. This statement is commented for the above reason.
 124 *
 125 *              INFR  =LINDEX(NLEVEL+1)
 126             ELSE
 127                INFR  =0
 128             ENDIF
 129             GO TO 100
 130          ELSE
 131 *
 132 *   *      Point is outside setup
 133 *
 134             NUMED = 0
 135             GO TO 999
 136          ENDIF
 137       ENDIF
 138 *
 139 *  **   Point is in current volume
 140 *
 141       IF(INFR  .GT.0) THEN
 142          JIN=LQ(JVO-INFR  )
 143          IQ(JIN) = IBSET(IQ(JIN),4)
 144          JVIN = JIN
 145       ENDIF
 146 * To avoid starting from the protuding part of a MANY volume
 147       IF(GONLY(NLEVEL).EQ.0.) THEN
 148         NLEVEL = NLEVEL -1
 149         GO TO 100
 150       ENDIF
 151       NLMIN = NLEVEL
 152       NLMANY = 0
 153 *
 154 * SECTION II: X is found inside current node at NLEVEL in /GCVOLU/.
 155 *             Search all contents recursively for any containing X.
 156 *             Take the first one found, if any, and continue at that
 157 *             level, incrementing NLEVEL and extending /GCVOLU/ tables.
 158 *             This is continued until a level is reached where X is not
 159 *             found in any of the contents, or there are no contents.
 160 * Note: Since Section II is re-entered from Section III, a blocking word
 161 * is used to mark those contents already checked.  Upon exit from Section
 162 * II, these blocking words are cleared at NLEVEL, but may remain set in
 163 * levels between NLEVEL-1 and NLMIN, if any.  They must be cleared at exit.
 164 *
 165 *  **  Check contents, if any
 166 *
 167   200 JVO = LQ(JVOLUM-LVOLUM(NLEVEL))
 168       NIN = Q(JVO+3)
 169       if(raytra.eq.1..and.imyse.eq.1)then
 170             CALL UHTOC(NAMES(NLEVEL),4,NAME,4)
 171             CALL GFIND(NAME,'SEEN',ISSEEN)
 172             if(isseen.eq.-2.or.isseen.eq.-1)goto 300
 173       endif
 174 *
 175 *   *   Case with no contents
 176 *
 177       IF (NIN.EQ.0) THEN
 178          GO TO 300
 179 *
 180 *   *   Case with contents defined by division
 181 *
 182       ELSEIF (NIN.LT.0) THEN
 183          CALL GMEDIV (JVO, IN, XC, 1)
 184          IF (IN.GT.0) THEN
 185             INFR   = 0
 186             GO TO 200
 187          ENDIF
 188 *
 189 *   *  Case with contents positioned
 190 *
 191       ELSE
 192        if(nin.gt.1)then
 193         clmoth=q(jvirt+4*(LVOLUM(NLEVEL)-1)+3)
 194         chmoth=q(jvirt+4*(LVOLUM(NLEVEL)-1)+4)
 195         ndivto=q(jvirt+4*(LVOLUM(NLEVEL)-1)+2)
 196         iaxis =q(jvirt+4*(LVOLUM(NLEVEL)-1)+1)
 197         if(iaxis.le.3)then
 198           ivdiv=((xc(iaxis)-clmoth)*ndivto/(chmoth-clmoth))+1
 199             if(ivdiv.lt.1)then
 200               ivdiv=1
 201             elseif(ivdiv.gt.ndivto)then
 202               ivdiv=ndivto
 203             endif
 204         else
 205           call gfcoor(xc,iaxis,cx)
 206           if(iaxis.eq.6)then
 207             if((cx-clmoth).lt.-1.)then
 208               cx=cx+360.
 209             elseif((cx-chmoth).gt.1.)then
 210               cx=cx-360.
 211             endif
 212             if(cx.gt.chmoth)then
 213               cx=chmoth
 214             elseif(cx.lt.clmoth)then
 215               cx=clmoth
 216             endif
 217           endif
 218           ivdiv=((cx-clmoth)*ndivto/(chmoth-clmoth))+1
 219             if(ivdiv.lt.1)then
 220               ivdiv=1
 221             elseif(ivdiv.gt.ndivto)then
 222               ivdiv=ndivto
 223             endif
 224         endif
 225         jvdiv=lq(jvirt-LVOLUM(NLEVEL))
 226         iofset=iq(jvdiv+ivdiv)
 227         ncont=iq(jvdiv+iofset+1)
 228         jcont=jvdiv+iofset+1
 229         if(ncont.eq.0)goto 260
 230        else
 231          JCONT  = LQ(JVO-NIN-1)+1
 232          NCONT  = 1
 233        endif
 234 *
 235 *     For each selected content in turn, check if point is inside
 236 *
 237          DO 259 ICONT=1,NCONT
 238            if(nin.eq.1)then
 239             in=1
 240            else
 241             IN = IQ(JCONT+ICONT)
 242            endif
 243             IF(IN.EQ.0) THEN
 244 *
 245 *     If the value IQ(JCONT+ICONT)=0 then we are back in the mother.
 246 *     So jump to 260, the search is finished. Clean-up should be done
 247 *     only up to ICONT-1, so we set:
 248 *
 249                NCONT=ICONT-1
 250                GOTO 260
 251             ELSE
 252             JIN = LQ(JVO-IN)
 253             IF (.NOT.BTEST(IQ(JIN),4)) THEN
 254                CALL GMEPOS (JVO, IN, XC, 1)
 255                IF (IN.GT.0) THEN
 256                   new2fl=0
 257                   IF (GONLY(NLEVEL).NE.0.) THEN
 258                     NLMANY = 0
 259                     nvmany = 0
 260                     nfmany = 0
 261                   ENDIF
 262                   INFR   = 0
 263                   GO TO 200
 264                ELSE
 265                   IQ(JIN) = IBSET(IQ(JIN),4)
 266                ENDIF
 267             ENDIF
 268             ENDIF
 269   259    CONTINUE
 270 *
 271   260    IF(NCONT.EQ.NIN) THEN
 272          DO 268 IN=1,NIN
 273             JIN = LQ(JVO-IN)
 274             IQ(JIN) = IBCLR(IQ(JIN),4)
 275   268    CONTINUE
 276          ELSE
 277          DO 269 ICONT=1,NCONT
 278            if(nin.eq.1)then
 279             in=1
 280            else
 281             IN  = IQ(JCONT+ICONT)
 282            endif
 283             JIN = LQ(JVO-IN)
 284             IQ(JIN) = IBCLR(IQ(JIN),4)
 285   269    CONTINUE
 286          IF(INFR  .GT.0) THEN
 287             JIN = LQ(JVO-INFR  )
 288             IQ(JIN) = IBCLR(IQ(JIN),4)
 289          ENDIF
 290          ENDIF
 291 *
 292       ENDIF
 293 *
 294 * SECTION III: X is found at current node (NLEVEL in /GCVOLU) but not in
 295 *              any of its contents, if any.  If this is a MANY volume,
 296 *              save it as a candidate best-choice, and continue the search
 297 *              by backing up the tree one node and proceed to Section II.
 298 *              If this is an ONLY volume, proceed to Section IV.
 299 *
 300 * *** Point is in current volume/medium, and not in any content
 301 *
 302   300 IF (GONLY(NLEVEL).EQ.0.) THEN
 303 *
 304 *  **   Lowest level is 'NOT ONLY'
 305 *
 306          IF (NLEVEL.GT.NLMANY) THEN
 307             CALL GSCVOL
 308             NLMANY = NLEVEL
 309             nfmany=nvmany+1
 310          ENDIF
 311          if(new2fl.eq.0)then
 312             nvmany=nvmany+1
 313             manyle(nvmany)=nlevel
 314             do 401 i = 1,nlevel
 315               manyna(nvmany,i)=names(i)
 316               manynu(nvmany,i)=number(i)
 317  401        continue
 318          endif
 319 *
 320 *   *   Go up the tree up to a volume with positioned contents
 321 *
 322          new2fl=-1
 323   310    INFR   = LINDEX(NLEVEL)
 324          NLEVEL = NLEVEL -1
 325          JVO    = LQ(JVOLUM-LVOLUM(NLEVEL))
 326          NIN    = Q(JVO+3)
 327          IF (NIN.LT.0) GO TO 310
 328 *
 329 C*****  Code Expanded From Routine:  GTRNSF
 330 C
 331          IF (GRMAT(10,NLEVEL) .EQ. 0.) THEN
 332             XC(1) = X(1) - GTRAN(1,NLEVEL)
 333             XC(2) = X(2) - GTRAN(2,NLEVEL)
 334             XC(3) = X(3) - GTRAN(3,NLEVEL)
 335 *
 336          ELSE
 337             XL1 = X(1) - GTRAN(1,NLEVEL)
 338             XL2 = X(2) - GTRAN(2,NLEVEL)
 339             XL3 = X(3) - GTRAN(3,NLEVEL)
 340             XC(1) = XL1*GRMAT(1,NLEVEL) + XL2*GRMAT(2,NLEVEL) + XL3*
 341      +      GRMAT(3,NLEVEL)
 342             XC(2) = XL1*GRMAT(4,NLEVEL) + XL2*GRMAT(5,NLEVEL) + XL3*
 343      +      GRMAT(6,NLEVEL)
 344             XC(3) = XL1*GRMAT(7,NLEVEL) + XL2*GRMAT(8,NLEVEL) + XL3*
 345      +      GRMAT(9,NLEVEL)
 346
 347          ENDIF
 348 C*****  End of Code Expanded From Routine:  GTRNSF
 349 *
 350          JIN = LQ(JVO-INFR  )
 351          IQ(JIN) = IBSET(IQ(JIN),4)
 352          NLMIN = MIN(NLEVEL,NLMIN)
 353          GO TO 200
 354       ENDIF
 355 *
 356 * SECTION IV: This is the end of the search.  The current node (NLEVEL
 357 *             in /GCVOLU/) is the lowest ONLY volume in which X is found.
 358 *             If X was also found in any of its contents, they are MANY
 359 *             volumes: the best-choice is the one among them at the greatest
 360 *             level in the tree, and it is stored.  Otherwise the current
 361 *             volume is the solution.  Before exit, all of the blocking
 362 *             words leftover in the tree must be reset to zero.
 363 * Note: A valid structure is assumed, in which no ONLY volumes overlap.
 364 * If this rule is violated, or if a daughter is not entirely contained
 365 * within the mother volume, the results are unpredictable.
 366 *
 367       DO 419 NL=NLMIN,NLEVEL-1
 368          JVO = LQ(JVOLUM-LVOLUM(NL))
 369          NIN = Q(JVO+3)
 370          DO 418 IN=1,NIN
 371             JIN = LQ(JVO-IN)
 372             IQ(JIN) = IBCLR(IQ(JIN),4)
 373   418    CONTINUE
 374   419 CONTINUE
 375 *
 376       if(nlmany.eq.0)then
 377         nvmany=0
 378         nfmany=0
 379       endif
 380       IF (NLMANY.GT.0) CALL GFCVOL
 381       JVO = LQ(JVOLUM-LVOLUM(NLEVEL))
 382       IF(JVIN.NE.0) IQ(JVIN) = IBCLR(IQ(JVIN),4)
 383       NUMED = Q(JVO+4)
 384 *                                                             END GMEDIA
 385   999 IF(JGSTAT.NE.0) CALL GFSTAT(2)
 386       END
 387 #endif